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FOREWORD 


This  document  is  the  Final  Report  on  research  on  “The  Numerical  Analysis  of  a  Class  of 
Problems  in  the  Mathematical  Theory  of  Plasticity  and  Damage.”  The  report  summarizes 
work  done  on  Contract  DAAL03-86-K-0043  between  the  Army  Research  Office,  Research 
Triangle  Park,  N.C.  and  The  University  of  Texas  at  Austin.  Principal  Investigator  of  the 
work  was  Professor  J.  Tinsley  Oden,  Director  of  the  Texas  Institute  for  Computational 
Mechanics.  This  report  provides  1)  a  list  of  publications  and  presentations  that  resulted  from 
the  work,  2)  a  list  of  students  and  faculty  supported  for  varying  periods  of  time  by  contract 
funds,  and  3)  chapters  summarizing  major  theoretical  and  numerical  results  obtained  during 
the  course  of  the  study. 
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SUMMARY  OF  WORK  AND  RESULTS 

1  Scope  and  General  Results  of  the  Project 

The  success  of  the  mathematical  sciences  in  providing  good  and  realistic  models  of  material 
behavior  has  led  to  great  expectations  in  the  engineering  community.  Today,  some  theories 
of  elasticity  and  plasticity  have  reached  a  level  of  maturity  that  admits  the  construction  of  a 
sound  and  rich  mathematical  basis  and  these  theories,  supported  and  confirmed  by  physical 
experiments,  in  turn  form  a  basis  for  the  numerical  analysis  of  many  problems  of  practical 
importance.  Other  theories  of  material  behavior  of  more  recent  origin  aim  at  characterizing 
more  subtle  or  more  complex  mechanical  phenomena  and, while  they  are  designed  with  the 
hope  of  modeling  many  fundamental  aspects  of  material  response,  they  are  often  founded  on 
shaky  mathematical  assumptions  which  block  further  development  and  implementation. 

Among  theories  of  material  behavior  that  fall  into  this  class  are  the  theories  of  damage 
that  grew  out  of  the  work  of  Kachanov  in  the  1950’s  and  which  have  fueled  a  growing 
percentage  of  the  engineering  literature  on  solid  mechanics  since  the  mid  1970’s.  Damage 
theory  purports  to  model  the  evolution  of  damage  of  a  material  body  subjected  to  forces  up 
to  the  initiation  of  macro  cracks  and  thus  to  provide  models  of  materials  just  prior  to  crack 
initiation  and  fracture.  There  does  not  exist  a  uniformly  accepted  definition  of  damage  in 
the  literature  (with  the  possible  exception  of  isotropic  damage)  and,  with  the  exception  of 
results  developed  in  the  present  study,  no  rigorous  mathematical  theory  of  damage  appears 
to  exist. 

One  characteristic  of  many  damage  theories  is  the  inclusion  of  an  evolution  equation 
for  the  growth  of  damage  in  a  material.  This  feature  establishes  a  formal  similarity  of 
damage  theory  to  modern  theories  of  viscoplasticity  and  rate-dependent  plasticity,  even 
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though  the  physical  mechanisms  responsible  for  damage  and  plasticity  are  quite  different. 
From  the  viewpoint  of  numerical  analysis,  many  numerical  techniques  effective  in  integrating 
the  ordinary  differential  equations  of  viscoplasticity  are  equally  effective  when  applied  to 
those  of  damage  theory.  Also,  as  was  shown  in  this  study,  techniques  can  be  developed 
which  are  applicable  to  problems  of  combined  viscoplastic  deformation  and  damage.  One 
result  of  numerical  experiments  done  in  this  study  was  the  discovery  that  many  (virtually 
all)  numerical  techniques  proposed  in  the  literature  for  handling  cyclic  viscoplasticity  and 
damage  are  unstable  or  only  marginally  stable,  and  special  steps  must  be  taken  to  produce 
reliable  numerical  simulations. 

It  was  these  issues  in  mind  that  the  present  investigation  was  initiated.  The  major  goals 
were: 

1.  the  study  and  formulation  of  a  physically  sound  and  mathematically  consistent  theory 
of  damage,  and,  in  particular,  to  explore  and  develop  the  qualitative  mathematical 
theory,  and 

2.  the  numerical  analysis  of  the  evolution  equations  of  damage  and  rate-dependent  plas¬ 
ticity. 

Some  new  developments  in  each  of  these  areas  were  contributed  in  the  project. 

2  General  Results 

A  fairly  detailed  study  of  many  existing  theories  of  damage  was  conducted.  One  conclusion, 
perhaps  not  surprising,  is  that  there  is  poor  agreement  among  reearchers  in  this  field  as 
to  what  consititutes  a  physically  correct  measure  of  damage  in  both  brittle  and  ductile 
materials.  For  anisotropic  damage,  scalar-,  vector-,  and  tensor- valued  damage  variables  have 
been  proposed  for  materials  that  undergo  elastic  and  elastoplastic  deformation.  Numerous 
deficiencies  and  inconsistencies,  both  physical  and  mathematical,  exist  in  some  of  the  more 
publicized  theories.  In  general,  the  field  is  still  quite  immature  and  the  general  acceptance 
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of  basic  principles  and  definitions  of  terms  have  neither  the  experimental  support  nor  the 
consensus  of  workers  in  the  field  to  form  the  nucleus  of  a  general  mathematical  theory. 

The  field  of  isotropic  damage  is  in  somewhat  better  shape.  In  the  present  study,  a  critical 
look  at  the  subject  was  conducted  and  several  new  results  were  produced.  In  addition  to 
providing  arguments  that  some  prevailing  theories  are  based  on  mathematically-unsound 
assumptions,  a  consistent  theory  for  isotropic  damage  was  identified  which  with  proper 
choices  of  parameters,  could  be  used  to  produce  results  consistent  with  experimental  data. 
An  interesting  and  possibly  new  result  of  this  study  was  the  establishment  of  actual  bounds 
on  the  moduli  and  Poisson’s  ratio  for  damaged  isotropic  elastic  solid.  Some  of  these  results 
are  summarized  in  Chapter  2  of  this  report. 

In  addition,  existence  and  uniqueness  theorems  were  proven  for  model  boundary-value 
problems  in  isotropic  damage  theory.  These  results  were  obtained  using  rather  standard 
techniques.  To  extend  them  to  more  general  theories  of  damage  and  plasticity  appears  to  be 
far  outside  the  reach  of  existing  mathematical  theories. 

A  second  part  of  the  study  focused  on  developing  numerical  models  of  damage  in  elastic- 
and  plastic  materials  and  in  developing  test  codes  for  numerical  experimentation.  In  this 
phase  of  the  work,  finite  element  models  of  damage  in  two-dimensional  structures  were 
developed  and  several  test  problems  were  solved.  The  most  challenging  feature  of  this  portion 
of  the  effort  was  the  development  of  robust  numerical  schemes  for  integrating  the  evolution 
equations  of  damage  and  rate-dependent  plasticity.  A  summary  of  the  work  on  this  subject  is 
reproduced  in  Chapter  3  of  this  report.  Briefly,  a  new  robust  implicit  scheme  was  developed 
which  proved  to  be  very  effective  in  treating  cyclic  problems  in  damage  simulations.  Studies 
on  error  estimation  and  stability  of  these  schemes  for  nonlinear  problems  remains  a  topic  for 
future  work. 
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3  Papers  and  Presentations 

Several  reports,  journal  articles,  and  presentations  resulted  from  work  on  this  project.  These 

are  listed  as  follows: 

1.  “A  Critique  and  Remarks  on  Damage  Theory,”  P.J.  Rabier  and  J.T.  Oden,  TICOM 
Report  82-10,  TICOM,  Austin,  1987. 

2.  “Some  Remarks  on  Damage  Theory,”  P.J.  Rabier,  International  Journal  of  Engineering 
Science, Volume  27,  No.  1,  pp.  29-54,  1989. 

3.  “Numerical  Solution  of  the  Evolution  Equations  of  Damage  and  Rate-Dependent  Plas¬ 
ticity,”  J.M.  Bass  and  J.T.  Oden,  International  Journal  of  Engineering  Science ,  Vol. 
26,  No.  7,  pp.  713-740,  1988. 

4.  “A  Note  on  the  Numerical  Analysis  of  Material  Damage  Based  on  the  Theory  of  Ma¬ 
terials  of  Type  N,”  S.J.  Kim  and  J.T.  Oden,  Computer  Methods  in  Mathematics  with 
Applications ,  Vol.  15,  No.  3,  pp.  169-174,  1988. 

5.  “Thermo- Viscoplastic  Analysis  of  Hypersonic  Structure  Subjected  to  Severe  Aerody¬ 
namic  Heating,”  Presented  at  the  30th  Structures,  Structural  Dynamics,  and  Materials 
Conference,  by  E.A.  Thornton,  J.T.  Oden,  W.  Tworzydlo,  and  S.K.  Youn,  Mobile,  AL, 
April  3-5,  1989. 

4  Students  and  Other  Personnel 

The  following  faculty  and  students  were  supported  for  varying  periods  of  time  on  the  project: 

1.  Dr.  J.T.  Oden,  Principal  Investigator 

2.  Dr.  P.J.  Rabier,  Visiting  Professor,  Mathematics,  University  of  Pittsburgh 

3.  Dr.  P.B.  Devloo,  Post-Doctoral  Research  Fellow 
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4.  Dr.  J.A.C.  Martins,  Visiting  Professor,  Technical  University  of  Portugal 

5.  Mr.  L.O.  Faria,  Ph.D.  student 


6.  Dr.  P.  Pattani,  Post-Doctoral  Research  Fellow 

The  following  students  received  some  support  from  the  project  for  a  short  period  of  time: 
Messrs.  0.  Hardy  and  T.  Westermann  and  Drs.  W.  Rachowicz,  C.Y.  Huang,  and  S.K.  Youn. 
Mr.  Hardy  completed  requirements  for  the  degree  of  Master  of  Science  and  L.O.  Faria  and 
W.  Rachowicz  completed  requirements  for  the  Ph.D.  degree  in  Engineering  Mechanics  during 
this  project.  While  their  theses  were  not  directly  related  to  the  theme  of  the  project,  they 
did  contribute  to  some  aspects  of  the  study  and  they  benefited  from  short-term  support  from 
the  contract. 
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SOME  REMARKS  ON  DAMAGE  THEORY 
Prologue 

This  Chapter  reproduces  the  studies  published  in  reports  and  papers  1,  2,  and  4 
listed  in  Section  1.2  of  the  previous  chapter.  It  is  intended  to  be  a  critique  of  several 
mathematical  aspects  of  damage  theories  proposed  in  contemporary  solid  mechanics 
literature. 

A  critical  look  at  this  literature  has  led  us  to  the  conclusion  that  to  date,  damage 
theory  suffers  from  both  the  lack  of  uniformity  in  its  approach  and  the  lack  of  rigor 
in  its  developments.  It  is  one  of  the  aims  of  this  paper  to  somewhat  remedy  the 
second  inconvenience  and  to  show  that,  in  turn,  rigor  greatly  helps  selecting  prevailing 
points  of  view.  This  programme  is  carried  out  for  the  broad  problem  of  anisotropic 
damage,  and  general  directions  for  future  research  in  various  areas  are  indicated.  In 
addition,  in  the  case  of  isotropic  damage,  we  have  been  able  to  establish  a  model  that 
does  not  exhibit  the  weaknesses  of  the  usual  generalization  of  Kachanov’s  theory. 
In  particular,  softening  of  the  material  is  translated  by  the  modifications  of  both 
Young’s  modulus  and  Poisson’s  ratio,  the  change  in  volume  due  to  voids  expansion  is 
properly  accounted  for,  and  the  need  for  empirically  determined  thresholds  for  failure 
is  eliminated.  The  model  is  valid  for  both  the  quasi-static  and  dynamic  cases,  in 
the  assumption  of  small  elasto-plastic  deformations  but  without  restriction  to  small 
damage.  For  simplicity,  attention  is  confined  to  isothermal  processes.  The  model  is 
well  suited  for  the  mathematical  treatment  of  P.D.E.’s.  The  final  section  is  indeed 
devoted  to  some  theoretical  existerce  and  uniqueness  results. 
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1  Introduction:  The  Origin  of  Damage  Theory 


Despite  that  the  macroscopic  notion  of  a  damaged  material  is  familiar  to  everyone’s 
most  primitive  intuition,  the  entire  process  through  which  an  originally  virgin  (un¬ 
damaged)  solid  becomes  damaged  is  far  from  transparent.  A  simple  but  useful  step 
towards  the  description  of  this  process  consists  in  making  explicit  the  common  under¬ 
standing  of  damage.  Usually,  damage  is  acknowledged  by  the  visual  observation  of 
cracks  which  were  not  present  in  the  virgin  material.  At  the  same  time,  the  relative 
size  and  the  number  of  such  cracks  will  automatically  be  taken  as  an  empirical  but 
instructive  estimate  on  the  level  of  damage. 

While  this  is  just  a  naive  statement  about  the  final  consequences  of  damage  which 
does  not  explain  at  all  how  such  cracks  have  been  able  to  develop,  nevertheless  it 
stresses  that  the  origin  of  damage  should  be  traced  back  to  properties  of  the  virgin 
material  that  are  capable  of  generating  such  cracks.  The  analysis  of  these  properties  is 
not  taken  up  by  fracture  mechanics,  devoted  to  the  study  of  evolution  of  pre-existing 
cracks,  but  which  does  not  investigate  the  structural  changes  occuring  between  the 
original  (virgin)  and  final  (cracked)  specimens. 

By  simply  ruling  out  the  existence  of  cracks,  classical  theories  such  as  elastic¬ 
ity  or  elasto-plasticity  equally  ignore  such  structural  changes.  The  reason  fracture 
mechanics  as  well  as  elasticity  or  elasto-plasticity  are  unable  to  provide  the  desired 
information  is  that  they  are  all  macroscopic  theories  whereas  damage  originates  in 
microphenomena  occuring  at  the  scale  of  grains.  These  phenomena  are  responsible 
for  local  modifications  of  the  elasto-plastic  properties  of  solids,  in  particular  elastic 
moduli,  during  1c  ling  and  unloading,  up  to  eventual  failure  (crack  initiation).  By 
establishing  p1:'>.'  les  appropriately  corroborated  by  experiments,  damage  theory, 
existing  or  to  be,  is  intended  to  provide  the  mathematical  models  of  this  transitional 
state. 

For  the  three  decades  since  the  pioneering  work  by  Kachanov  [14],  the  starting 
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point  of  damage  theory  has  been  experimental  evidence  that  any  parcel  of  material 
contains  a  multitude  of  defects  in  the  form  of  microvoids.  During  a  loading- unloading 
process,  these  voids  may  undergo  irreversible  growth  and  new  ones  may  nucleate.  The 
ultimate  coalescence  of  these  voids  results  in  macro-crack  initiation.  With  Krajcinovic 
[17],  we  note  that  nucleation  and  growth  of  such  microvoids  is  a  mechanism  distinct 
from  propagation  of  dislocations  responsible  for  plasticity.  Just  as  in  the  macroscopic 
case,  the  number  and  size  of  voids  in  any  unstressed  configuration  may  be  taken  as  an 
indication  on  the  current  level  of  damage.  The  crucial  difference  is  that  the  smaller 
scale  now  under  consideration  enables  one  to  use  these  measurable  data  as  parameters 
representing  damage  at  any  given  material  point.  Rather  unfortunately  however,  no 
decisive  argument  has  been  made  yet  that  dictates  a  specific  way  to  measure  damage. 
One  of  the  many  problems  is  the  vague  reference  to  some  unstressed  configuration 
that  generally  does  not  exit.  This  will  be  clarified  below  and,  on  the  basis  of  a  short 
critical  review,  we  shall  see  that  a  point  of  view  definitely  seems  to  prevail. 

2  The  Measure  of  Damage 

In  what  follows,  “damaged  material”  refers  to  a  certain  specimen  having  undergone 
a  certain  loading/unloading  process  during  a  specified  period  of  time.  By  “virgin 
material”,  we  mean  the  same  specimen  before  any  such  process  has  started. 

In  his  original  one-dimensional  model,  Kachanov  (loc.  cit.)  suggests  that  damage 
should  be  defined  as  the  density  of  microvoids  in  any  cross  section  of  the  material. 
More  specifically,  consider  a  rod  submitted  to  uniform  tractions  and/or  compressions 
in  the  direction  of  its  axis,  so  that  the  deformation  of  the  rod  can  accurately  be 
described  by  a  one-dimensional  model.  A  mathematical  point  X  is  then  associated 
with  any  cross  section  of  the  rod.  The  relative  area  of  the  voids  in  the  cross  section 
at  X  defines  a  scalar  variable  D(X)  (0  <  D(X )  <  1)  representing  damage. 

In  the  three-dimensional  case,  one  can  attempt  to  construct  a  simple  generalization 
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of  this  idea  as  follows:  let  X  be  a  point  of  the  three-dimensional  material  in  its 
current  configuration  and  consider  a  small  volume  element  AV  about  X.  Since  the 
geometry  of  the  microvoids  in  AV  is  altered  by  the  elastic  part  of  the  deformation, 
AF  must  first  be  put  in  a  stress-free  configuration,  which  can  be  done  by  separating 
AV  from  the  remainder  of  the  body  and  freed  of  external  forces.1  We  shall  henceforth 
denote  by  AV'  the  unstressed  configuration  of  AF,  unique  to  within  rigid  motions, 
and  assume  after  suitable  translation  that  the  transformation  AV  — *  AF  leaves  X 
invariant  (this  is  done  for  notational  convenience  only).  Each  plane  through  X  with 
normal  n  intersects  AV  along  a  planar  surface  A/4,  depending  on  n.  The  relative 
area  of  microvoids  in  AA,  say  D(X ,  n)  can  be  taken  as  a  measure  of  damage.  If  the 
distribution  of  voids  is  isotropic,  D(X,  n)  is  independent  of  n  and  becomes  a  variable 
D(X)  as  in  Kachanov’s  model.  Mathematically,  considering  the  limiting  case  of  an 
infinitesimal  volume  element  dF,  the  above  procedure  amounts  to  defining  D{X)  as 
a  density  of  microvoids. 

Although  this  approach  has  been  used  even  in  recent  work  in  the  literature,  this 
choice  for  the  damage  variable  suffers  several  criticisms.  First,  although  if  it  is  rea¬ 
sonable  to  assume  that  the  distribution  of  voids  is  isotropic  in  the  virgin  material, 
there  is  clear  experimental  evidence  that  this  isotropy  is  destroyed  as  a  result  of  pla¬ 
nar  voids  (microcracks)  nucleating  and  growing  “in  the  direction  perpendicular  to 
the  maximum  tensile  strain”  (Krajcinovic  and  Fonseka  [18],  referring  to  experimental 
work  by  Hayhurst  and  Leckie).  Isotropic  damage  can  therefore  be  expected  only  as 
a  result  of  isotropic  stress  fields,  a  situation  of  little  practical  importance.  But  there 
are  a  few  exceptions  to  this  rule  and  some  materials,  in  particular  some  aluminum 
alloys,  do  exhibit  essentially  isotropic  damage  more  or  less  independently  of  the  stress 
field.  In  such  cases,  it  is  also  clear  from  experimental  observations  that  damage  takes 
the  form  of  nucleation  and  growth  of  nearly  spherical  microcavities.  Because  it  ncc- 

*At  the  scale  of  the  volume  AV,  it  is  reasonable  to  assume  zero  stress  in  absence  of  external 
forces;  this  follows  from  the  postulate  that  material  points  retain  no  residual  stress,  unlike  the  whole 
body  if  plastic  effects  have  taken  place. 
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essarily  goes  along  with  modification  of  volume  through  the  principle  of  conservation 
of  mass  (unlike  damage  due  to  microcracks)  this  phenomenon  is  better  described  in 
terms  of  the  relative  volume  of  the  microcavities.  As  we  shall  later  see,  no  purely 
mathematical  argument  can  be  made  to  link  the  two  methods  of  measurement. 

A  third  criticism  is  that,  in  the  framework  of  the  theory  of  small  elasto- plastic 
deformations,  the  damage  variable  D(X )  defined  above  is  used  jointly  with  the  notion 
of  effective  stress  (stress  relative  to  the  effective  area  1  —  D{X))  to  yield  a  questionable 
model  in  which  softening  of  the  material  occurs  without  modification  of  Poisson’s 
ratio. 

In  accordance  with  the  foregoing  arguments  against  the  previous  definition  for 
the  damage  variable,  one  may  arrive  at  the  conclusion  that  the  geometrical  difference 
between  planar  voids  (microcracks)  and  microcavities  is  reflected  by  two  different 
aspects  of  damage:  on  the  one  hand,  nucleation  and  growth  of  microcavities  will  be 
accompanied  by  changes  in  volume  (and  mass  density)  consistent  with  macroscopic 
measurements  (porosity)  and  will  not  generate  anisotropy  in  the  material.  Anisotropy 
will  be  accounted  for  by  nucleation  and  growth  of  planar  voids  which,  in  turn,  will  not 
induce  modification  of  volume.  It  is  therefore  both  natural  and  desirable  to  represent 
these  two  aspects  of  damage  by  two  complementary  variables. 

The  obvious  choice  for  defining  a  damage  variable  representing  nucleation  and 
growth  of  microcavities  is  their  volume  density,  denoted  by  fl(.X’)  at  each  point  X 
of  the  material.  For  practical  measurements,  this  density  can  be  identified  with  the 
relative  volume  of  microcavities  in  the  previous  volume  element  AK  so  that,  again, 
0  <  f2(X)  <  1.  For  a  virgin  material,  £l(X)  =  0,2  i.e.,  is  small  enough  to  be 
negligible.  As  a  result,  if  x  denotes  the  position  originally  occupied  by  X  in  the 
virgin  material  and  pa(x )  is  the  mass  density  of  the  virgin  material  at  x,  then  the 

2This  assumption  is  not  essential  and  the  initial  condition  can  be  taken  as  an  arbitrary  function 
with  values  in  [0, 1].  If  so,  formula  (2.1)  must  be  modified  accordingly. 
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mass  density  p(X )  of  the  damaged  material  at  X  would  be  given  by 

p(X)  =  (1  -  fl(X))p.(x),  (2.1) 

if  the  current  configuration  were  stress-free. 

The  above  method  of  measuring  isotropic  damage  due  to  microcavities  through 
the  scalar  variable  fi  was  introduced  and  used  by  Davison,  Stevens  and  Kipp  [8],  and 
one  can  hardly  imagine  a  drastically  different  and  nonequivalent  procedure.  But  the 
approach  for  measuring  damage  due  to  microcracks  is  far  from  being  as  self-evident. 

The  question  here  is  to  account  for  both  the  area  of  the  cracks  and  their  orien¬ 
tation.  The  importance  of  the  orientation  is  obvious  from  the  intention  of  relating 
microcracks  to  anisotropy.  It  has  often  been  taken  as  a  corollary  to  the  direction 
dependence  that  the  corresponding  damage  variable  should  not  be  a  scalar:  a  vector 
representation  has  been  proposed  by  Davison  and  Stevens  [7]  in  which  the  direction 
of  the  vector  is  an  average  of  the  directions  of  the  normals  to  the  cracks  lying  in  the 
volume  element  AV,  and  its  magnitude  is  the  projected  area  of  the  cracks  onto  the 
plane  orthogonal  to  the  average  direction. 

This  definition  is  ambiguous,  for  it  is  not  an  easy  task  to  average  directions  (these 
being  equally  represented  by  two  opposite  unit  vectors).  In  particular,  the  only 
consistent  way  to  define  the  average  direction  of  cracks  uniformally  distributed  is  to 
choose  it  to  be  the  zero  vector  (!)  and  there  are  other  problems,  such  as  obtaining  a 
reasonably  smooth  vector  field  as  a  result  of  the  operation  of  averaging  directions.  To 
partly  circumvent  these  difficulties,  the  authors  quoted  above  have  suggested  the  use 
of  not  one  but  several  vector  fields,  each  one  of  them  accounting  for  cracks  with  more 
or  less  similar  orientation  in  the  undamaged  material.  Aside  from  its  rather  empirical 
character,  this  approach  prohibits  taking  the  original  damage  to  be  zero  and,  more 
generally,  does  not  allow  for  a  correct  evaluation  of  damage  caused  by  nucleation  of 
new  cracks  if  their  direction  is  not  already  included  in  the  list  of  vector  fields. 

Serious  concern  should  be  expressed  on  the  basis  that  nucleation  appears  to  be  a 
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dominant  factor  during  most  of  the  damaging  process  (Onat  and  Leckie  [25])  whereas, 
as  pointed  out  before,  orientation  of  cracks  depends  on  the  applied  forces.  Another 
objection  is  made  by  Krajcinovic  et  al.  [17,18]  who  observed  that  the  model  does  not 
reduce  to  Kachanov’s  for  uniaxial  problems,  the  difference  being  the  area  of  the  planar 
voids  entering  the  constitutive  equations  linearly  (in  [14])  instead  of  quadratically  (in 
[8]).  While  an  ad-hoc  adjustment  is  proposed  in  [17,18]  to  reconcile  both  points  of 
view,  the  very  problem  of  definition  of  the  vector-valued  damage  variable  is  prudently 
eluded  by  making  rapid  reference  to  “an  appropriate  averaging  process.  .  .  ”  with  no 
further  detail. 

Undoubtedly,  the  surviving  advocates  of  vector  representations  are  only  motivated 
by  the  desire  of  simplifying  other  existing  tensor  representations.  Regarding  these, 
we  shall  only  mention  that  there  are  several  points  of  view  expressed  in  the  literature, 
as  to  what  the  tensor’s  coefficients  should  represent  and  what  its  order  should  be. 
Not  all  of  them  are  free  of  criticism  either.  The  interested  reader  may  find  relevant 
references  and  additional  comments  in  the  aforementioned  papers  by  Krajcinovic  et 
al.,  but  it  is  our  feeling  that  some  significant  progress  has  recently  been  made  by 
Onat  and  Leckie  [25].  Their  approach,  with  a  technical  difference  in  the  method  of 
measuring  damage  (cf.  Remark  2.2  later)  is  as  follows. 

To  begin  with,  let  us  move  one  step  backwards  and  deny  the  existence  of  planar 
voids.  This  merely  means  that,  no  matter  how  flat  it  may  look,  a  void  has  a  positive 
volume  or  else  is  not  a  void.  We  are  then  ready  to  repeat  the  procedure  leading  to 
the  density  D(X,  n)  introduced  earlier.  This  variable  is  still  scalar  valued  but,  for 
fixed  X ,  is  a  function  of  the  direction  n.  Equivalently,  D(X,  •)  identifies  with  an 
(even)  function  on  the  unit  sphere  52  and,  as  such,  can  be  expressed  as  a  generalized 
Fourier  series  by  the  means  of  spherical  functions.  Taking  advantage  of  the  evenness 
of  D(X,  •)  and  after  a  suitable  grouping  of  terms,  one  arrives  at 

D(X,n)  =  D0(X)+j2kHDi](X)+  £  fijU(n)  Dijkl(X)  +  •  •  • ,  (2.2) 


namely,  an  expression  involving  even  rank  tensors  of  basis  functions  (1,  (/,j(n)), 
(fijki(n)),  •••)  together  with  even  rank  tensors  of  coefficients  (1, (Dij(X)), 
(Dijki(X)),  . . . ).  The  tensors  of  basis  functions  have  explicit  expressions  in  terms 
of  the  direction  n;  for  instance,  in  obvious  notation 

ftj (n)  =  TliTlj  1  ^  <  3, 

and  the  tensors  of  coefficients  can  be  obtained  through  integral  formulae  (in  which 
integration  is  performed  w.r.t.  n  €  52),  e.g. 

D0(X)  =  i/ '  D(X,n), 

JS2 

OiiiX )  =  g  D(X,n)fij(n),  l<ij<3.  (2.3) 

From  (2.3)  it  is  clear  that  D0(X)  represents  the  isotropic  contribution  since 
D0  =  D  whenever  D(X,  n)  is  independent  of  n  while  all  other  coefficients  Dij{X ), 
Dijki{X),  ...,  vanish  due  to  uniqueness  of  the  decomposition  (2.2).  All  the  terms 
Djj(X),  Dijki{X),  in  (2.2)  therefore  translate  the  anisotropy  of  the  distribu¬ 
tion  of  voids.  Of  course,  for  practical  purposes,  one  needs  to  assume  that  only  a  few 
terms  are  significant  in  (2.2)  throughout  the  loading/unloading  process:  starting  with 
isotropic  damage  (e.g.,  zero)  in  the  virgin  material,  one  may  assume  that 

D(X,n)  =  D0(X), 

for  materials  complying  with  the  isotropic  damage  assumption,  the  simplest  case  of 
anisotropic  damage  being  handled  through  the  choice 

D(X,  n)  =  D0(X)  +  £  fij{n)Dij(X).  (2.4) 

0=1 

Because  of  symmetry  (Dij  =  Dji)  and  traceless  property  (Dn  +  D2 2  +  D33  =  0) 
(cf.  [15])  the  above  formula  shows  that  the  simplest  model  for  anisotropic  damage 
involves  six  scalar  variables  instead  of  three  in  a  vector  representation.  Whether  or 
not  acceptable  accuracy  is  achieved  with  (2.4)  certainly  depends  on  the  material  and 
the  applied  forces.  If  not,  one  knows  immediately  from  (2.2)  how  to  include  terms 
providing  a  better  approximation. 
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Remark  2.1:  The  above  considerations  do  not  justify  the  choice  of  the  expansion 
(2.2)  versus  the  usual  series  of  spherical  functions.  The  reason  has  to  be  found  in  the 
fact  that  each  term  in  (2.2)  involves  an  even  rank  tensor,  say  2 p  with  p  =  0,  1,  ... , 
of  coefficients.  This  tensor  must  be  viewed  as  being  p  times  covariant  and  p  times 
contravariant.  As  such,  we  are  prompted  to  introduce  a  state  variable  since  the  way 
this  term  complies  with  changes  of  coordinates  is  compatible  with  a  state  variable 
being  independent  of  the  coordinate  frame.  For  more  detail,  see  Fardshisheh  and 
Onat  [9].  By  comparison,  the  series  of  spherical  functions  equivalent  to  (2.2)  follows 
from  the  decomposition  of  the  space  of  functions  on  S?  into  subspaces  invariant  under 
the  action  of  the  group  of  rotation  (see  Gel’fand  and  Shapiro  [10]).  As  opposed  to 
any  of  their  particular  bases,  these  spaces  are  the  important  intrinsic  entities  that 
capture  anisotropy.  However,  they  only  yield  a  decomposition  of  D(X,n)  as  a  series 
of  vectors  which  are  not  appropriate  state  variables. 

Remark  2.2:  In  [25],  where  the  above  approach  is  introduced,  the  choice  for  the 
damage  variable  is  different.  Its  definition  is  based  on  the  optical  observation  that 
in  materials  such  as  copper  exhibiting  anisotropic  damage,  microvoids  nucleate  and 
grow  on  grain  boundaries.  Assuming  the  grain  boundaries  to  be  planar  and  given 
a  unit  vector  n  (E  S2,  consider  a  surface  element  AS  C  S2  centered  at  n.  Picking 
X  and  AV  as  before,  call  D(X,n )  the  total  volume  of  voids  on  grain  boundaries 
orthogonal  to  some  n!  €  A 5,  relative  to  AS  (i.e.,  divided  by  the  area  of  AS)  in  the 
stress-free  configuration  AV  of  AV.  In  the  limiting  case  of  an  infinitesimal  surface 
element  dS,  this  defines  D(X,n )  as  a  density  of  voids  volume  on  grain  boundaries 
per  unit  area  of  S2 .  However,  a  slight  modification  is  needed  if  one  wants  to  preserve 
the  desirable  property  that  the  quantities  D0,  Dij ,  etc.  ...in  (2.3)  are  densities 
themselves.  Indeed,  in  the  present  formulation,  one  has  /  D(X,  n)  ~  2Vj  where 
Vj  is  the  total  volume  of  voids  in  AV.  Hence,  Da(X )  =  Vyj'ln  is  not  a  density, 
and  difficulties  arise  when  considering  the  limiting  case  of  an  infinitesimal  volume 
element  dV .  But  this  is  easily  overcome:  it  suffices  to  replace  D(X,n)  above  by 
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D(X  ,n)/\AV\  where  jAl^l  denotes  the  volume  of  AV  (so  that  D(X,n)  appears  as 
a  density  per  unit  area  of  S 2  of  a  density  of  voids  per  unit  volume!).  In  particular, 
D0(X)  becomes  D0(X)  =  V^-/2n|AV|,  indeed  a  density,  and  the  same  is  true  of  the 

other  coefficients  Dij,  etc . This  definition  turns  out  to  have  also  an  advantage 

regarding  kinematical  considerations  (see  Section  4).  Relating  anisotropic  damage  to 
voids  on  grain  boundaries  leaves  no  other  option  than  relating  isotropic  damage  to 
voids  nucleating  and  growing  within  the  grain.  This  attitude  seems  to  be  in  perfect 
agreement  with  theoretical  and  experimental  results.  For  completeness,  it  must  be 
mentioned  that  some  materials  (e.g.,  steel)  exhibit  both  kinds  of  void  distributions. 


When  damage  is  isotropic,  it  is  tempting  to  try  to  relate  D(X,n)  =  D(X)(= 
D0(X)  in  (2.3))  to  the  other  variable  fI(J£)  introduced  for  the  same  purpose.  That 
this  is  not  possible  through  purely  mathematical  arguments  can  be  seen  as  follows: 
pick  AV  such  that  AV  is  a  spherical  ball  centered  at  X  and  consider  any  plane 
through  X  intersecting  AV  along  the  disk  AA.  From  the  definitions  of  D  and  f l,  one 
has,  modulo  higher  order  terms, 


D(X)  = 


A 

|£*f 


where  A  and  V  denote  the  area  and  volume  of  voids  in  AA  and  A  V,  respectively,  and 
where  |A/l|  and  \AV\  stand  for  the  area  and  volume  of  AA  and  AV,  respectively. 
Clearly 

9tt  D3{X)  A 3 

16  W(X)  ~  Ve  ' 

Thus,  if  fl(X)  can  be  expressed  as  a  function  of  D(X),  one  finds,  for  a  given  value  of 
D(X),  a  formula  for  the  calculation  of  A  in  terms  of  V.  To  see  that  such  a  formula 
does  not  exist,  consider  two  holes  in  AV  having  volume  V,  but  one  occupying  a 
spherical  ball  centered  at  X  and  the  other  a  spherical  layer  along  the  boundary  of 
AV.  In  both  cases,  A  is  independent  of  the  plane  through  X  but  A  =  |A4|fl(X)2/3  in 
the  first  case  while  A  =  |A<4|(1  —(1  — fl(A"))2^3)  in  the  second  case.  The  impossibility 
of  obtaining  a  formula  for  A  in  terms  of  V  (or  conversely)  is  easily  traced  back  to  the 
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fact  that  euclidian  measure  in  spherical  coordinates  depends  on  the  distance  to  the 
origin,  and  hence  the  position  and  geometry  of  the  holes  cannot  be  left  unspecified 
in  the  calculation  of  their  total  volume  from  the  knowledge  of  A .  Also,  it  is  not 
hard  to  convince  oneself  that  the  counter-example  given  above  represents  the  two 
extreme  cases  of  any  distribution  of  voids  in  AV.  It  follows  that  D(X)  can  always 
be  estimated  from 

1  -  (1  -  n(X))7/3  <  D(X)  <  Q(X)2/3.  (2.5) 

3  Damage  and  Elastic  Properties 

We  shall  begin  with  some  notation:  let  U  denote  the  original  (t  —  0)  unstressed 
configuration  of  the  virgin  material,  identified  with  an  open  subset  of  R3  and  let  <f)(U) 
be  the  configuration  at  some  time  t  =  t0  >  0.  In  the  configuration  (f)(U ),  the  body  is 
in  a  certain  state  of  stress  and,  during  the  time  interval  [0,  t0\,  has  undergone  loading 
and  unloading  in  a  way  that  may  have  caused  inelastic,  permanent  deformations. 

Through  nucleation  and  growth  of  microvoids,  damage  is  responsible  for  at  least 
part  of  these  inelastic  deformations.  Therefore,  unless  damage  has  not  evolved  and  no 
other  plastic  phenomenon  has  taken  place  in  the  time  interval  [0,  i„],  the  deformation 
4>  is  not  an  elastic  one.  In  standard  plasticity  theory,  one  would  assume  that  the 
elastic  properties  of  the  material  in  the  configuration  4>{U)  are  the  same  as  those  of 
the  material  in  configuration  U.  This  means  that  the  constitutive  equation  between 
stress  and  elastic  strain  is  unaffected  by  the  elastic  deformation.  In  damage  theory, 
the  bottom  line  is  that  the  elastic  properties  of  the  material  have  been  modified  as 
a  result  of  damage.  A  major  component  of  the  study  is  then  the  establishment  of 
sound  models  for  these  modifications. 

In  the  assumption  that  damage  is  a  state  of  the  material,3  damage  theory  natu¬ 
rally  fits  into  the  framework  of  elasto-plasticity  with  state  variables  developed  since 
3Some  theories  do  not  start  with  this  assumption;  see  Section  4. 
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the  mid-sixties  for  various  purposes.  Apparently,  nothing  has  been  done  in  the  always 
difficult  case  of  finite  deformations.  Even  the  more  tractable  problem  of  infinitesimal 
deformations  is  quite  unsatisfactorily  solved,  mainly  because  the  abundance  of  defini¬ 
tions  for  the  damage  variable  makes  impossible  joint  efforts  in  a  uniformly  accepted 
direction. 

Although  it  is  apparently  always  taken  for  granted,  we  believe  that  the  preliminary 
step  of  investigating  how  the  modifications  of  the  elastic  properties  can  be  measured  is 
worth  a  few  comments.  Certainly,  the  elastic  properties  of  the  material  at  a  point  X  of 
the  configuration  <f>(U )  cannot  be  measured  directly  from  this  configuration.  Neither 
can  they  be  measured  from  a  configuration  obtained  from  <t>(U)  after  unloading  since, 
even  if  damage  and  plastic  deformation  are  supposed  to  be  unaffected  by  the  unloading 
process,  the  unloaded  configuration  will  generally  not  be  free  of  residual  stress.  As 
was  the  case  with  the  measure  of  damage,  one  must  take  advantage  of  the  fact  that 
the  elastic  properties  of  the  body  at  X  are  not  dependent  on  the  entire  configuration 
<f>(U)  but  can  (and  must)  be  determined  after  an  elementary  volume  A  V  about  X  has 
been  separated  from  the  remainder  of  the  configuration  <£(£/)  and  freed  of  external 
forces  so  that  it  is  now  stress  free  in  the  configuration  AV.  Again,  we  have  used  the 
postulate  that  material  points  retain  no  residual  stress. 

The  deformation  between  AV  and  AV  is  an  elastic  one,  at  least  under  the  rea¬ 
sonable  assumption  that  unloading  occurs  with  no  change  in  plastic  deformation  (nor 
damage,  but  this  is  obvious  from  the  very  definition  for  the  measure  of  damage  made 
in  Section  2).  Assuming  also  this  deformation  to  be  infinitesimal,  the  stress  at  X  in 
the  configuration  AV  (i.e.,  4>(U))  is  a  linear  function  of  the  strain  tensor  of  AV  —>  AV 
at  X  (the  configuration  AV  being  chosen,  as  before,  so  as  to  leave  X  invariant).  Ob¬ 
serve  in  passing  that  in  the  limiting  case  of  an  infinitesimal  volume  dV ,  this  strain 
tensor  coincides  with  the  so-called  elastic  strain  of  the  deformation  <f>  at  X.  As  is 
well  known,  the  determination  of  the  elastic  response  of  AV  depends  on  21  con¬ 
stants.  In  this  respect,  it  is  an  implicit  hypothesis  of  damage  theory  that  anisotropy 
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(resp.  isotropy)  of  the  distribution  of  voids  in  AV  is  responsible  for  anisotropy  (resp. 
isotropy)  of  the  elastic  response.  The  one  and  only  difficulty  is  then  to  make  the 
relationship  explicit. 

Even  the  isotropic  case  in  which  the  problem  reduces  to  finding  the  two  Lame 
constants  is  not  agreed  upon.  For  easy  reference,  we  shall  speak  of  “Lame  constants 
at  X ”  although  the  above  arguments  clearly  show  that  they  are  not  calculated  in 
the  configuration  4>(U).  The  model  for  isotropic  damage  based  on  an  extension  of 
Kachanov’s  approach  discussed  in  the  previous  section  assumes  that  both  Lame  con¬ 
stants  at  X  =  <f)(x),  x  £  U,  are  obtained  by  multiplication  of  the  Lame  constants  at 
x  4  through  the  same  factor  1  —  D(X).  The  invariance  of  Poisson’s  ratio  mentioned 
earlier  follows  immediately.  The  models  based  on  vector  representation  of  the  dam¬ 
age  variable  lead  to  a  five-constant  constitutive  equation  (assuming  there  is  only  one 
damage  field)  reminiscent  of  those  employed  in  the  assumption  of  transverse  isotropy. 
In  these  models,  the  Lame  constants  are  unaffected  by  damage  which  is  only  involved 
in  the  other  three.  The  rigorous  method  used  by  Onat  and  Leckie  to  measure  damage 
has  still  to  be  accompanied  by  a  thorough  examination  of  how  the  various  (tensor) 
variables  should  enter  the  constitutive  equation. 

To  be  complete,  let  us  mention  that  thermal  effects  cannot  be  ignored  in  a  some¬ 
what  general  theory  of  damage,  although  their  study  may  reasonably  be  postponed 
until  the  isothermal  case  is  better  understood.  Thermal  effects  are  actually  included 
in  the  work  of  Davison  et  al.  [8]  on  isotropic  damage,  in  which  an  interesting  idea 
of  how  to  calculate  the  Lame  constants  at  X  is  introduced.  In  the  assumption  of 
negligible  thermal  effects,  we  shall  now  expand  on  a  generalization  of  this  idea  and 
derive  some  interesting  and  unexpected  consequences. 

The  starting  point  is  the  simple  observation  that  both  the  Lame  constants  at  X 
and  the  damage  variable  D(-X’)  introduced  in  Section  2  are  calculated  from  the  same 
'’Here,  the  terminology  applies  in  a  strict  sense  since  U  is  supposed  to  be  a  natural  state 
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volume  element  AV  about  A”  in  a  stress-free  configuration  AV.  The  microvoids  in 
AV  thus  occupy  the  volume  fraction  ft(A)  and  hence  the  volume  fraction  1  —  0(A) 
only  consists  of  actual  material.  If  now  AV  is  considered  as  a  composite  microvoids- 
material,  there  are  several  available  procedures  to  estimate  the  Lame  constants  of 
the  composite.  As  in  [8],  we  shall  follow  the  self-consistent  method  of  Budiansky  [3] 
which  yields,  since  the  Lame  constants  of  a  microvoid  must  obviously  be  taken  equal 
to  zero 


K 

V 


K, 


Vo 


3(1  —  v) 

2(1  —  2u) 

15(1  -  u) 


n 


7  -5v 


n 


(3.1) 

(3.2) 


where  K  and  K0  (resp.  fi  and  fi0)  are  the  bulk  moduli  (resp.  shear  moduli)  at  A  and 
x  =  A)  respectively,  u  =  u(X)  is  Poisson’s  ratio  at  A  (i.e.,  of  the  composite) 
and  ft  stands  for  ft(A).  Naturally,  v  is  also  related  to  K  and  /x  through 


3  K  -  2  n 
2(3  A'  +  /O’ 


(3.3) 


a  formula  complementing  (3.1)  and  (3.2)  and  allowing  for  actual  calculation  of  K 
and  /x  in  terms  of  ft.  Note  that  (3.1  )-(3.3)  only  yield  implicit  characterization  of  K 
and  //;  note  also  that  the  conditions  /x  >  0  and  A  >  0  (with  A'  =  A  +  |/x)  are  not 
guaranteed  by  the  above  relations  and  depend  on  ft  which  then  should  range  in  an 
interval  [0,  ftmajt)  with  ftmax  <  1. 


Although  its  analog  for  arbitrary  composite  materials  cannot  be  solved  explicitly 
in  general,  the  system  (3.1  )-(3.3)  does  admit  an  explicit  solution.  This  can  be  seen 
as  follows:  setting 

0  =  K  In,  (3.4) 


one  finds 


30  —  2 
2(30+1)' 


(3.5) 
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For  future  use,  it  is  also  worth  noticing  that 


A  >  0,  n  >  0  <=>  0  >  -,  p  >  0 


From  (3.5),  we  get 


'  1  -2u  =  3/(30  +  1), 

•  1  -i/  =  (30  +  4)/2(30  +  1), 

_  7-5z/  =  (270  +  24)/2(30+l). 

With  60  =  K„/Ho ,  it  follows  from  (3.5)  and  (3.7)  that 

(4  —  30f2  —  4O)(90  +  8) 
%(90  +  8  -  150Q  -20fl)  ' 


(3.6) 


(3.7) 


(3.8) 


Hence,  0  can  be  recovered  from  0O  and  fl  by  solving  a  quadratic  equation,  namely 


(36  -  60H  +  270oCl)02  +  (32  -  80H  -  360O  +  600ofl)0  +  320o(tt  -  1 )  =  0  (3.9) 

It  is  easily  checked  that  this  equation  has  exactly  one  positive  root  for  Q  (E  [0, 1]  if 
0O  >  8/9. 

Let  g{ fl)  denote  the  positive  root  of  (3.9).  Substituting  0  —  g(fl)  into  (3.7)  and 
from  (3.1)  and  (3.2),  one  finds  the  desired  expression  for  K  and  p  as  functions  of 
fi.  Of  course,  they  are  too  complicated  to  be  of  any  use  in  theoretical  arguments. 
However,  to  corroborate  the  intuitive  idea  that  the  material  becomes  softer  and  softer 
as  the  damage  Q  grows,  it  is  natural  to  expect  I\  and  g  to  be  decreasing  functions  of 
fi.  Also,  it  would  be  of  interest  to  determine  the  value  flmax  beyond  which  either  A 
or  g  becomes  negative,  thus  ceasing  to  represent  an  actual  material.  These  questions 
can,  in  fact,  be  answered  with  practically  no  calculation  if  one  uses  (3.8)  to  derive 
the  expression  of  0  in  terms  of  0  and  0o,  which  is  easily  seen  to  be 

St  -  no t  =  1  (g.-Q)(M  +  8) 

n  '~3  (0  +  4/3)((90„  —  20)0  +  80„]  1  ’ 

A  key  observation  is  that,  irrespective  of  0o  ^  4/3,  not  only  the  substitution  0  =  0o 
in  (3.10)  yields  fl  =  0  but  the  substitution  0  =  4/3  yields  Q  =  1/2.  In  addition,  it  is 


21 


immediate  to  check  that  for  9  =  4/3  and  fl  =  1/2,  both  constants  K  and  g  vanish. 
This  suggests  that  the  study  of  Q  =  f(Q)  as  given  by  (3.10),  but  restricted  to  the 
interval  /  =  [0o,4/3j  if  0o  <  4/3,  /  =  [4/3,^]  if  90  >  4/3,  will  provide  the  desired 
information. 

To  begin  with,  note  that  /  is  well  defined  on  I  since  ( 9a  ^  4/3  is  always  implicit 
in  what  follows) 

80o/(90o  -  20)  <  90  if  90  <  4/3, 

80o/{99o  -  20)  >  60  if  4/3  <90<  20/9, 

8 90/{990  -  20)  <  0  if  90  >  20/9, 

and  there  is  no  problem  if  90  —  20/9.  It  is  equally  straightforward  to  check  that 
/(/)  C  [0, 1/2]  irrespective  of  60 ,  hence 

/(/)  =  [0, 1/2], 

from  f(0o)  =  0,  / (4/3)  =  1/2  and  the  intermediate  value  theorem.  Moreover,  for 
every  9  £  I  and  fl  =  /(<?),  9  is  uniquely  characterized  to  be  the  unique  positive 
root  g( Cl)  of  equation  (3.9).  This  shows  that  /  and  g  are  the  inverse  of  each  other. 
Actually,  this  has  been  shown  for  90  >  8/9  only,  the  only  case  for  which  equation  (3.9) 
has  exactly  one  positive  root  for  every  fl  €  [0, 1].  But  since  the  interval  of  interest  is 
[0, 1/2]  instead  of  [0, 1]  as  would  have  been  expected  at  the  beginning,  and  since  it  is 
straightforward  to  check  that  equation  (3.9)  has  a  unique  positive  root  irrespective 
of  fl  G  [0, 1/2]  for  every  90  >  0,  the  restriction  90  >  8/9  is  not  necessary.  From  the 
fact  that  /  and  g  are  inverse  of  each  other,  it  follows  that  /  is  a  bijection  from  /  to 
[0, 1/2].  Continuity  being  obvious,  we  infer  that  /  is  monotone  on  /  and,  finally,  that 
g  =  f~l  is  monotone  on  [0, 1/2]  s. 

Suppose  first  that  g  is  increasing.  Using  (3.1)  and  (3.7),  we  see  that 

I<{=  A'(fi))  =  7v0(l  -  l(g(Q)  +  l)ft),  (3.11) 

which  immediately  shows  that  I\  is  a  decreasing  function  of  0  G  [0, 1  / 2].  In  particular 
K( f!)  >  0  for  fi  G  [0, 1/2)  since  A'(l/2)  =  0  is  already  known. 

5Trying  to  give  a  direct  proof  of  monotonicity  of  /  and  g  for  every  ^4/ 3  is  rather  inextricable. 
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Next,  from  (3.4)  we  find 


0  =  =  A'(n)/5(ft). 

Since  K  >  0  and  p  >  0  on  [0, 1/2],  K  is  decreasing  and  (?  is  increasing,  it  follows  that 
g  is  decreasing.  In  particular,  g(fi)  >  0  for  Cl  6  [0,1/2)  since  //( 1/2)  =  0  is  already 
known.  From  the  equivalence  (3.5)  and  recalling  that  0o  >  2/3  in  physically  relevant 
situations,  hence  <?(0)  >  min(0o,4/3)  >  2/3  for  every  0  g  [0, 1/2],  we  obtain 

Omax  =  1/2, 


irrespective  of  0o  ^  4/3. 

Suppose  now  that  g  is  decreasing.  Using  (3.2)  and  (3.7),  we  see  that 

20/3 


/*(= 


\  / 5  ,  20/3  ^ 

-/i°l  \3  +9ff(n)  +  sJ 


a 


(3.12) 


95(0)  +  8; 

which  immediately  shows  that  g  is  a  decreasing  function  of  0  E  [0, 1/2].  In  particular, 
/<(0)  >  0  for  ft  £  [0,1/2)  since  /z(l/2)  =  0  is  already  known.  Next,  from  (3.4)  we  find 


I<(=  I<(Q))  =  g(Sl)g(  0). 

Since  g  >  0  and'5  >  0  on  [0, 1/2]  and  g  and  g  are  both  decreasing,  it  follows  that  K 
is  decreasing.  Assuming  again  6a  >  2/3,  the  conclusion 


Omax  =  1/2, 

is  reached  through  the  same  arguments  as  above. 

For  0  =  1/2,  we  know  that  the  corresponding  value  of  9  is  0  =  4/3.  With  (3.5), 
we  deduce  the  corresponding  value  of  Poisson’s  ratio 

1 

V  ~  5' 

The  above  analysis  is  valid  under  the  assumption  0o  ^  4/3.  For  0o  =  4/3  (i.e. 
v0  =  1/5  if  v0  denotes  Poisson’s  ratio  of  the  virgin  material),  equation  (3.9)  has 
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9  =  4/3  as  its  unique  positive  root  regardless  of  ft  €  [0,1]  (so  that  the  relation 
between  6  and  ft  is  not  invertible  and  the  previous  approach  does  not  go  through). 
Substituting  9  =  4/3  into  (3.7),  relations  (3.1)  and  (3.2)  become 

K  =  A'0(l  —  2ft),  /i  =  /i0(l-2  ft), 

and,  again,  we  find 

Ana*  =  1/2. 

Because  9  =  4/3  is  independent  of  ft,  no  modification  of  Poisson’s  ratio  (=  1/5)  is 
observed  in  this  limiting  case. 

Remark  3.1:  More  generally,  as  a  result  of  (3.5),  it  follows  that  v  =  i/(ft)  is 
decreasing  (resp.  increasing)  when  g  is  decreasing  (resp.  increasing).  As  v(  1/2)  =  1/5 
irrespective  of  initial  conditions,  it  follows  that  g  is  decreasing  if  v0  >  1/5  (i.e., 
90  >  4/3),  increasing  if  uQ  <  1/5  (i.e.,  90  <  4/3). 

The  value  ftmaK  =  1/2  instead  of  the  expected  value  1  is,  after  all,  not  surprising. 
Let  us,  for  instance,  decide  to  view  <f>(U)  as  a  union  of  identical  microcubes,  each  of 
them  containing  the  same  spherical  hole.  It  is  obvious  that  the  structure  will  be  unable 
to  support  any  stress  when  the  volume  of  the  spherical  hole  is  maximum,  namely 
when  it  occupies  7r/6(~  0.52)  of  the  volume  of  the  cube.  The  relation  ftmax  =  0.5  can 
roughly  be  explained  on  this  basis. 

Let  us  point  out  that  considering  the  damaged  material  as  a  composite  to  which 
the  rules  dictated  by  self-consistent  methods  apply  is  not  the  only  possible  attitude. 
The  proposed  model  may  require  modifications  for  values  of  damage  approaching 
1/2,  based  either  on  refinements  of  the  theory  or  empirical  observations.  But  the 
fact  that  it  is  self-sufficient  to  predict  vanishing  of  the  Lame  constants  before  the 
maximum  possible  value  1  for  the  damage  variable  is  important.  This  phenomenon, 
interpreted  as  failure  (crack  initiation)  is  indeed  well  known  but  is  not  included  in 
simpler  models  which  must  appeal  to  experimentally  determined  thresholds.  In  [20] 
Lemaitre  indicates  the  values  0.2  to  0.8  for  the  range  of  critical  values  relative  to  the 
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variable  D(X),  depending  on  the  material.  Existence  of  such  thresholds  is  not  in 
contradiction  with  the  relation  Qmax  =  1/2  regardless  of  the  material.  Indeed,  recall 
that  there  is  no  mathematical  correlation  between  D(X)  and  fl(.X’),  except  for  the 
estimate  (2.5).  In  practice,  this  means  that  two  materials  showing  isotropic  damage 
may  perfectly  exhibit  different  maximum  values  D^x  of  the  variable  D  despite  that 
their  values  of  fimax  is  the  same.  This  is  because  the  average  size  of  each  individual 
void  may  be  typically  different  in  the  two  materials  for  specimens  of  equal  total 
volume. 

As  a  matter  of  fact,  it  is  not  even  certain  that  Dma *  is  independent  of  the  ex¬ 
periment  for  a  given  material.  Substituting  f2(X)  =  fimax  =  1/2  into  (2.5)  yields 
the  bracket  [0.37,0.63]  for  Z)max  This  bracket  is  narrower  than  [0.2, 0.8]  previously 
reported.  However,  it  is  not  clear  that  the  latter  has  actually  been  established  with 
materials  that  all  show  isotropic  damage,  which  could  explain  the  discrepancy. 

Earlier  in  this  section,  we  have  established  that  both  A'(f2)  and  fi(fl)  are  decreasing 
functions  of  ft.  This  was  done  to  corroborate  the  fact  that  the  material  was  indeed 
softening  a s  damage  increases.  Later,  in  connection  with  the  dissipative  aspect  of 
damage  evolution,  we  shall  need  the  following  less  intuitive  and  stronger  result: 

Proposition:  Set 

it  (a)  =  (i  -si)-'/3k(Q), 
m)  =  (i  —  n)-,/V(o)- 

Then ,  K(Cl)  and  p(U)  are  descreasing  functions  of  fl  G  [0, 1/2]. 

Proof:  Again,  we  shall  distinguish  the  two  cases  when  the  function  g(Vt)  giving 

the  unique  positive  root  of  equation  (3.9)  is  either  increasing  or  decreasing  on  [0, 1  /2] 
(recall  that  there  is  no  other  option). 

Suppose  then  that  g  is  increasing.  From  (3.11)  and  a  straightforward  calculation, 
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one  finds 


/<•'(«)  <  o  <=>  -f(a'(n)n  +  s(fi)  +  i)(i  -  n>  +  |[i  -  |(9(n)  +  i)n]  <  o. 

The  second  inequality  can  be  rewritten  as 

3!'  - ?(s(«)  +  OT  <  j(s'(Si)n  +  9(fi)  +  i)(i  -  n), 

and,  since  g'(£l)  >  r  will  certainly  be  satisfied  if 

|[i  -  +  i)n]  <  n)  +  i)(i  -  n)  <=> 

3  <  (g(Sl)  +  i)(f  - 

As  fl  G  [0, 1/2],  this  inequality  follows  from 

3  <  I(5(0)  +  1), 

which  is  obvious  because  g  >  0.  Hence,  A'(fi)  is  a  decreasing  function  of  H. 

To  show  that  /i(fi)  is  decreasing  too,  let  us  simply  observe  that  A'(fl)//i(fi)  = 
K(Cl)/n(Sl)  =  g(£l),  so  that  /x(H)  =  A'(fl)/<?(H).  The  result  is  immediate  from  K 
being  decreasing  and  nonnegative  and  g  being  increasing  and  nonnegative. 

Suppose  now  that  g  is  decreasing.  With  (3.12),  one  finds 

m)  <  0  <=>  N\U)(  1  -Q)  +  IJV(Q)  <  0,  (3.13) 
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The  second  inequality  in  (3.13)  is  then  true  provided  that 

N(Q)  <5(1  -fl). 

But  this  is  immediate  for  it  is  obvious  that  N(Cl)  <  1  while  5(1  —  ft)  >  5/2  for  Cl  £ 
[0,1/2].  Hence,  (i(Cl)  is  a  decreasing  function  of  Cl  Writing  as  before  K  (Cl)  /  fi(Cl)  = 
g(Cl),  namely  K(Cl)  =  ft(Cl)g(Cl),  one  finds  that  K(Cl)  is  decreasing  as  the  product  of 
two  decreasing  nonnegative  functions  and  the  proof  is  complete. 

4  Kinematics 

From  now  on,  x  will  denote  an  arbitrary  point  of  the  reference  configuration  U  and 
X  =  <t>(x)  represents  the  same  material  point  in  the  current  configuration.  This  is 
consistent  with  notation  used  in  previous  sections.  If  the  problem  were  considered 
in  the  framework  of  elasto-plasticity  theory,  the  deformation  gradient  V<^(*)  would 
be  decomposed  according  to  one  of  several  available  procedures  into  “elastic”  and 
“plastic”  parts.  Most  commonly,  one  would  use  the  multiplicative  decomposition  of 
Lee  [19]  to  write 

V0(*)  =  E(x)P(x),  (4.1) 

where  E  and  P  denote  elastic  and  plastic  components,  respectively.  Because  dislo¬ 
cations  merely  introduce  negligible  variations  of  volume,  it  is  required  that 

det  P(x)  =  1.  (4.2) 

The  decomposition  (4.1)  thus  does  not  allow  for  the  inelastic  dilatation  typical  of 
isotropic  damage.  Incidentally,  note  that  this  observation  alone  justifies  that  damage 
be  considered  a  phenomenon  distinct  from  plasticity. 

In  [8],  Davison  et  al.  have  made  use  of  a  modification  of  the  decomposition  (4.1) 
first  introduced  by  Kratochvil  for  other  (but  somewhat  related)  purposes,  and  written 
V<£(sc)  in  the  form 

V</>(jb)  =  E(x)M(x)P(x),  (4.3) 
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where  now  the  tensor  M{x)  accounts  for  the  effect  of  damage.  It  is  fairly  simple  to 
determine  the  structure  of  the  tensor  M(x).  For  this,  it  is  helpful  to  refer  once  again 
to  the  volume  element  AV  about  X  =  <}>(x)  and  its  stress  free  configuration  AV  (see 
Section  3).  Clearly,  considering  the  limiting  case  of  an  infinitesimal  volume  element 
dV,  one  sees  that  E(x )  coincides  with  the  gradient  of  the  (elastic)  deformation  be¬ 
tween  dV  and  dV.  To  determine  M(x),  consider  the  volume  element  Av  =  <f>~x(AV) 
about  x.  Except  for  a  plastic  deformation  which  does  not  affect  its  volume,  AV  dif¬ 
fers  from  Av  through  the  dilatation  introduced  by  evolution  of  damage.  Since  At;  and 
AV  contain  the  same  volume  of  material,  which  represents  only  the  fraction  1  —  fl(X) 
in  AV,  one  has  (assuming  zero  damage  in  the  reference  configuration)  modulo  higher 
order  terms 


|AV|  = 


|Au| 


i -si{xy 

Considering  the  limiting  case  of  infinitesimal  elements,  we  find  that,  necessarily 


det  M(x)  =  (1  —  Q(J£))-1,  X  =  <f>(x).  (4.4) 


As  the  distribution  of  microvoids  is  isotropic  and  M{x)  depends  only  on  this  distribu¬ 
tion,  no  direction  must  be  privileged  in  M(x)  (i.e.,  the  first  order  difference  between 
Av  and  AV  consists  of  a  stretch  equal  in  all  directions).  It  follows  that  M(x )  is  a 
multiple  of  identity  and,  together  with  (4.4),  we  find 

M(x)  =  (1  -  fl(X))-1/3/,  X  =  <f>(x)  .  (4.5) 

From  now  on,  it  will  be  more  convenient  to  introduce  the  damage  variable 

u>(x)  =  f }(<t>(x )),  X  6  U,  (4.6) 

expressed  in  the  reference  configuration.  Combining  (4.3)  and  (4.5),  we  arrive  at 

V<j>{x)  =  (1  -  u(x))-l/3E(x)P(x).  (4.7) 
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Setting 

<f>(x)  =  X  +  u(x), 

<  E(x)  =  I  +  e(x),  (4.8) 

.  P(x)  =  I  +  p(a J), 

an  equivalent  form  of  (4.7)  is 

Vit(x)  =  [(1  -  u{x))~^3  -  1] I  +  (1  -  a>(a:))-1/3(e(a:)  +  p{x)  +  e(x)p(x)). 

In  the  case  of  infinitesimal  elastic  and  plastic  deformations,  this  relation  reduces  to 
(omitting  x  -  dependence) 

Vu  =  [(1  -  u>)~1/3  -  1]/  +  (1  -  u;)~1/3(e  +  p). 

Therefore,  using  the  standard  notations 

e(u )  =  |(Vm  +  VuT), 

•  ee  =  |(e  +  eT)> 

.  ep  =  |(P  +  pT), 

we  find 

e(u)  =  [(1  -  u)-1'3  -  1] I  +  (1  -  w)-1/3(ee  +  ep) 

or,  equivalently 

ee  =  (1  -  w)1/3e(«)  -  [1  -  (1  -  u)l'3}I  -  ep  . 

Relation  (4.2)  is  taken  into  account  by  requiring 

Trsp  =  0.  (4.12) 

If  anisotropic  damage  were  to  be  considered,  the  decomposition  (4.3)  would  still 
be  appropriate,  but  the  tensor  M(x)  should  also  represent  the  anisotropy  of  the 
distribution  of  voids,  hence  could  no  longer  be  a  multiple  of  identity  in  general. 
However,  repeating  previous  arguments,  its  determinant  should  still  relate  to  the 
total  volume  of  voids.  This  demonstrates  a  significant  advantage  in  the  choice  of  the 


(4.9) 

(4.10) 

(4.11) 
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density  of  voids  volume  on  grain  boundaries  per  unit  area  of  the  unit  sphere  as  a 
measure  of  damage  (cf.  Remark  2.2).  Indeed,  recall  that  the  corresponding  quantity 
Da(X)  in  (2.2)  relates  to  the  total  volume  Vp  of  voids  in  AK  through  the  simple 
relation 

D0{X)  =  Vp/2n\AV\  . 

Hence 

det  Af  (*)  =  (1  -2wD0(X))~1. 

If  the  damage  variable  D(X,  n)  is  defined  as  in  Section  2,  an  experimentally  deter¬ 
mined  correlation  between  Da(X)  and  Vp  must  be  used,  assuming  of  course  that 
such  a  correlation  exists  for  a  given  material.  In  this  respect,  recall  that  no  purely 
mathematical  argument  can  be  made  to  establish  the  desired  correlation.  In  any 
case,  M(x )  will  depend  on  the  damage  variables  accounting  for  the  anisotropy  of  the 
distribution  of  voids.  The  difficulty  in  relating  M(x)  to  these  variables  is  greatly 
dependent  on  the  clarity  of  their  geometrical  interpretation  (but  the  problem  is  nev¬ 
ertheless  not  expected  to  be  trivial).  In  this  respect  too,  the  variables  introduced 
through  the  approach  of  Onat  and  Leckie  seem  to  be  the  most  appropriate. 

Remark  4.1:  In  the  event  that  anisotropy  in  the  distribution  of  voids  can  be 
attributed  to  the  plastic  slip  alone,  the  tensor  M(x)  remains  a  multiple  of  identity. 
Such  an  assumption  may  only  apply  to  materials  for  which  damage  is  isotropic  when 
it  occurs  in  the  elastic  range.  We  have  found  no  information  in  the  literature  as  to 
whether  these  materials  form  a  broader  class  than  those  for  which  damage  is  always 
isotropic.  For  the  elements  of  this  class  and  for  them  only,  anisotropy  is  entirely 
translated  by  the  nature  of  the  elastic  response  (constitutive  equation).  For  the 
record,  we  note  that  it  has  generally  been  overlooked  that  anisotropy  must  also  be 
involved  in  the  kinematical  aspect  of  the  problem. 

To  complete  this  section,  we  shall  briefly  mention  how  the  decomposition  (4.3) 
can  be  used  as  a  starting  point  for  a  totally  different  approach  to  damage  theory. 
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Let  us  begin  with  the  assumption  that,  due  to  damage,  the  tensor  V<fr(x)  has  the 
decomposition  (4.3)  instead  of  (4.1).  For  a  virgin  material,  M(x)  =  J.  This  suggests 
writing  M(x)  =  I  +  m(x)  and  hence,  omitting  se-dependence. 

/  +  2 s(u)  +  VmtVu  =  I  +  2ee  +  2em  -f  2ep  +  higher  order  terms,  (4.13) 

where 

£m  =  ^(m  +  mr), 

appears  as  a  damage  strain.  Assuming  infinitesimal  strains,  relation  (4.13)  becomes 

e(u)  =  e‘  +  em  +  ep.  (4.14) 

Decomposition  (4.14)  above  is  used  e.g.  by  Nicholson  [24],  Evolution  of  damage  is 
next  translated  through  a  suitable  differential  equation  for  the  damage  strain  em,  just 
as  plastic  phenomena  require  the  establishment  of  an  appropriate  differential  equation 
for  the  plastic  strain  ep.  This  approach  does  not  require  a  specific  damage  variable 
to  be  introduced,  but  has  the  concomitant  disadvantage  of  not  relating  damage  to  its 
physical  origin  (on  the  positive  side  is  the  fact  that  one  does  not  have  to  deal  with 
anisotropy,  as  a  result  of  the  voids  not  being  explicitly  involved 

in  the  formulation  of  the  problem).  Such  a  theory  however  can  hardly  be  used 
to  investigate  the  modifications  of  the  elastic  properties  due  to  damage.  These  mod¬ 
ifications  must  then  be  negligible  for  the  validity  of  the  theory.  In  turn,  negligible 
modification  of  the  elastic  properties  enforces  the  hypothesis  of  small  damage. 

5  A  Model  for  Isotropic  Damage 

The  loading  and  unloading  of  the  initially  virgin  material  will  be  taken  into  account  by 
time-dependent  body  force  densities  6  (per  unit  mass)  and  surface  traction  densities 
t  (per  unit  area)  defined  in  the  current  configuration  <f>((J).  During  this  process, 
supposed  to  be  isothermal ,  the  body  is  clamped  along  the  part  T0  C  T  =  dU  of  the 
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reference  configuration  U.  In  particular,  the  surface  tractions  t  are  prescribed  on 
0(I\)  where  Ti  =  r\ro.  The  total  body  force  acting  at  X  €  4>{U)  is  the  sum 

Div  ^(T^)  +  pb, 

where  T^(=  T^{X))  is  the  Cauchy  stress  tensor  at  X ,  Div^  is  the  divergence 
operator  in  the  coordinate  frame  of  the  deformed  configuration  and  p  denotes  mass 
density  at  X.  For  X  G  <£(I\),  the  total  surface  traction  vanishes,  i.e., 

(T^  +  t)dA^  =  0, 

with  and  dA &  being  the  unit  outward  normal  vector  and  the  infinitesimal  area 
element  at  X,  respectively.  Using  the  Piola  transformation,  i.e.,  setting 

T(x)  =  detV0(x)T^(«£(x))V0(x)-T,  x  =  <t>~\X),  (5.1) 

allows  one  to  derive  corresponding  expressions  in  the  reference  configuration  U 
through  the  identity  (cf.  Ciarlet  [4]). 

Div  T(x)  =  det  V<t>{x)  Div  ^(T^X^x)),  x  G  U, 

T(x)v(x)dA  =  T^>(</>(x))i/<^(0(a;))dA^>, 

where  Div  is  the  divergence  operator  in  the  coordinate  frame  of  the  reference  con¬ 
figuration  and  i/(x)  and  dA  are  the  unit  outward  normal  vector  and  infinitesimal 
area  element  at  x,  respectively.  Thus,  the  total  body  force  at  X  =  <f>(x)  is  (using 
po(x)  =  p(X)  det  V<£(x)) 

[det  V<£(x)]-1[  Div  T(x)  +  p0(x)b(<f>(x))],  x  G  U,  (5.2) 

while  balance  of  surface  tractions  reads 

-  T(x)v(x)dA  =  t{<j>{x))dA $  =  0,  x  G  ra.  (5.3) 

In  the  hypothesis  of  dead  loading ,  the  body  force  density  b(<f>(x))  and  the  surface 
traction  densities  are  independent  of  <j>,  i.e., 

b(4>(x))  =  6(x),  t(^>(x))d/l^>  =  t(x)dA. 
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In  this  case,  it  easily  follows  from  (5.2)  and  (5.3)  that  balance  of  linear  momentum 
yields  (omitting  ^-dependence) 

Div  T  -I-  pab  =  p0u  in  t/, 

-  u  =  0  on  r„,  (5.4) 

Tv  —  t  on  Tj, 

where  we  have  used  <f>  =  I  +  u. 

We  shall  limit  ourselves  to  considering  the  case  of  infinitesimal  elasto-plastic  de¬ 
formations.  In  the  notation  of  Section  3,  the  elastic  part  of  the  deformation  can 
be  viewed  as  the  deformation  between  the  volume  elements  AV  and  AV.  For  an 
infinitesimal  elastic  deformation,  this  means  that 

T&  =  A(u)ee ,  (5.5) 

where  u>  is  defined  in  (4.6)  and,  since  the  natural  configuration  AU  is  isotropic  and 
isothermal  conditions  are  assumed,  the  fourth-order  tensor  A(u)  is  defined  by 

A(u>)e  =  A(u >)(Tre)I  +  2p(uj)e,  (5.6) 

for  every  3x3  (symmetric)  tensor  e.  Of  course,  the  Lame  constants  A(u>)  and  p(u>) 
will  be  taken  according  to  the  procedure  of  Section  3. 

Using  (5.1)  together  with  the  hypothesis  of  infinitesimal  elasto- plastic  deforma¬ 
tions  in  (4.7)  and  (4.8),  the  system  (5.4)  may  then  be  rewritten  as  (neglecting  higher 
order  terms) 

f  Div  [(1  -  u;)~2/3A(u;)ee]  +  pQb  =  p„u  in  17, 

<  u  =  0  on  r0, 

(  (1  —  u>)-2/3.4(u;)ee  —t  on  Tj  . 

In  turn,  recalling  (4.11),  this  becomes 

Div  [B(w)(£(u)  -  (1  -  u>)-1/3ep)]  —  V7r(u;)  +  pQb  =  p0it  in  U, 

-  u  =  0  on  ro,  (5.7) 

[j3(cj)(e(w)  —  (1  —  u>)-1/,3eP)j  v  =  7r(u>)»/  -f  t  on  I\, 

where 

B(u)  =  (1  -  w)-1/3>l(u>),  (5.8) 
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and  where  the  pressure  tc  depends  only  on  damage  through  the  relation 

t(u»)  =  3A'(w)(l  -  w)',/3[(l  -  u)-1'3  -  1],  (5.9) 

with  K(u)  =  A(w)  +  In  particular,  ir(u)  vanishes  at  failure  (ui  =  1/2), 

since  K{u))  does  (cf.  Section  3).  This  should  be  expected  from  the  fact  that  x(u>)I  is 
essentially  the  residual  stress  (Remark  5.2). 

Remark  5.1:  It  is  important  to  observe  that  despite  the  assumption  of  infinitesimal 
elasto-plastic  deformations,  the  strain  tensor  e{u)  need  not  be  small.  This  has  indeed 
not  been  required  in  the  foregoing  analysis  and  such  an  assumption  would  contradict 
relation  (4.10)  since  the  coefficient  of  /  can  be  as  large  as  0.26  ( for  u  =  1  /2).  Besides, 
on  the  basis  of  (4.10)  again,  requiring  e(u)  to  be  small  would  immediately  limit 
damage  to  be  small,  too.  This  shows  the  importance  of  not  having  such  a  restriction. 

Remark  5.2:  Suppose  that  motion  has  occurred  only  in  the  elastic  regime  (i.e., 
ep  =  0)  and  the  forces  have  been  constant  long  enough  for  the  body  to  be  in  equilib¬ 
rium  in  the  configuration  </>(£/).  Then,  the  system  (5.8)  reduces  to 

—  Div  (B(w)e(u))  =  b  —  V7r(u>)  in  U, 

<  u  =  0  on  ro, 

(B(u)e(u))i/  =  7r(u>)i/  +  t  on  Tj. 

This  system  is  consistent  with  a  model  of  linear  elasticity  in  which  the  reference 
configuration  U  is  isotropic  but  is  not  a  natural  state  under  zero  external  forces,  unless 
7 r(w)  =  0,  i.e.,  lj  =  0  or  1/2  at  each  point.  This  agrees  with  the  physical  argument 
that  void  nucleation  and  growth  demands  an  increase  of  volume  for  a  natural  state 
(if  any  exists  at  all)  and  the  body  cannot  be  confined  to  the  volume  available  in  U 
without  the  action  of  a  residual  stress.  Consistency  with  isotropy  follows  from  the  fact 
this  this  residual  stress  is  here  the  pressure  n(u>)I.  Naturally,  except  when  ro  =  T, 
the  reference  configuration  U  is  not  an  actual  configuration  of  the  material  in  general. 

For  our  model  to  be  complete,  we  must  introduce  suitable  evolution  equations  for 
the  plastic  strain  ep  and  for  the  damage  variable  u.  A  differential  equation  of  the 
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form 

w  =  F^(u;,r^),  (5.10) 

where  F &  is  a  function  to  be  determined  from  experiments,  is  of  common  use  in 
the  literature.  This  is  what  is  done  e.g.,  by  Trampczynski,  Hayhurst  and  Leckie 
[27],  Davison  et  al.  [8].  The  former  work  uses  a  damage  variable  u  defined  from 
D  (cf.  Section  2)  instead  of  fi,  namely  u:(x)  =  £>(<£(  sc)).6  But  one  must  be  aware 
that  a  formulation  like  (5.10)  cannot  be  valid  without  severe  limitations,  because  it 
is  not  compatible  with  the  basic  postulate  that  damage  is  not  affected  by  (at  least 
appropriate)  unloading.  Recall  that  this  postulate  enters  the  very  definition  of  the 
damage  variable  (cf.  Section  2).  Rather,  one  should  say  that  (5.10)  holds  when 
loading  occurs,  while  u>  =  0  during  unloading.  The  ambiguity  is  of  course  to  define 
loading  versus  unloading,  which  can  be  done  only  in  the  case  of  proportional  loading. 
Therefore,  it  appears  that  a  relation  such  as  (5.10)  can  be  valid  only  for  monotone 
proportional  (or  nearly  proportional,  although  “monotone”  must  then  be  understood 
in  some  empirical  sense)  loading.  This  limited  framework  is  nevertheless  relevant  in 
many  practical  situations. 

In  [8],  where  no  distinction  among  stresses  needs  to  be  made  in  view  of  the  com¬ 
pletely  linear  model  adopted  (in  particular,  ignoring  the  residual  stress),  the  function 
F&  in  (5.10)  is  derived  from  a  model  for  nucleation  and  growth  of  microvoids  devel¬ 
oped  by  Barbee,  Seaman,  Crewdson  and  Curran  [2].  In  this  respect,  let  us  point  out 
that  it  is  indeed  worthwhile  taking  the  hint  from  the  materials  science  approach  to 
derive  an  appropriate  formula  for  the  rate  u.  However,  individual  void  growth  is  often 
measured  through  parameters  which  do  not  directly  relate  to  the  density  u ;  (see  e.g. 
Cocks  and  Ashby  [6])  and  some  difficulties  arise  in  establishing  an  appropriate  cor¬ 
relation.  Also,  an  important  ingredient  in  the  determination  of  F $  which  is  omitted 
by  studies  on  individual  void  growth  is  coalescence  of  two  or  more  voids:  the  rapid 
yielding  of  the  ligament  between  two  adjacent  voids  resulting  in  coalescence  will  make 
6As  has  been  the  case  throughout  this  paper,  time-dependence  is  implicit. 


35 


the  damage  rate  u>  larger  than  what  it  would  be  in  a  situation  when  voids  grow  in 
a  relatively  independent  fashion,  at  least  when  coalescence  is  significant.  Obviously, 
the  possibility  of  coalescence  increases  with  the  number  of  voids.  This  agrees  with 
the  premise  that  nucleation  is  dominant  at  the  beginning  of  damage  while  growth 
takes  over  during  later  stages,  as  noted  by  Onat  and  Leckie  [25]. 

Remark  5.3:  Coalescence  of  two  voids  is  often  referred  to  as  “rupture”  or  “failure” 
by  materials  science  specialists.  This  definition  is  somewhat  questionable  since  coa¬ 
lescence  of  two  voids  merely  provides  a  bigger  void,  whose  size  is  however  at  the  same 
scale  as  the  size  of  the  original  two.  In  our  definition,  failure  occurs  when  the  voids 
density  u>  is  large  enough  to  make  the  elastic  moduli  vanish  (u>  =  1/2  if  the  approach 
of  Section  3  is  taken).  This  makes  failure  relate  to  a  property  of  the  material,  with 
no  concern  as  to  how  the  critical  value  of  u>  is  reached  (nucleation,  individual  growth, 
coalescence  or  any  combination  of  these). 

Summarizing  the  preceding  comments,  the  function  F&  should  then  be  taken  as 
the  sum  of  the  three  rates 

nucleation  -f  individual  growth  +  coalescence.  (5-11) 

For  the  first  term,  one  may  refer  to  the  work  by  Goods  and  Brown  [11]  (see  also 
Argon,  Im  and  Safoglu  [1]).  Individual  growth  from  the  point  of  view  of  materials 
science  has  already  been  mentioned  above.  Both  the  first  two  terms  are  considered 
in  Barbee  et  al.  (loc.  cit.)  who  do  treat  void  growth  through  the  volume  density 
u>.  But  their  theory  is  admittedly  limited  to  small  damage  and  underestimates  the 
actual  void  volume  by  a  factor  of  2  for  larger  values  of  w.  We  strongly  suspect  that 
this  is  due  to  not  considering  the  third  term  in  (5.11),  for  which  further  investigation 
is  needed.  It  should  be  negligible  for  small  damage  and  dominant  as  u>  approaches 
1/2. 

We  observe  that  an  expression  that  agrees  with  (5.11)  has  been  used  by  Tvergaard 
[28],  although  not  in  a  damage  model  in  the  sense  of  this  paper,  and  according  to 
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other  approaches  (see  the  references  in  (28]).  These  approaches  lead  to  a  function 
F&  in  (5.10)  which  also  depends  on  the  stress  rate  T^\  As  we  shall  later  see,  such 
expressions  are  useful  to  deal  with  nonproportional  or  cyclic  loading.  Allowing  F & 
in  (5.10)  to  depend  on  ep,  possible  dependence  on  the  plastic  strain  rate  kp  can  also 
be  incorporated  in  the  hypothesis  of  rate-independent  plasticity,  namely  assuming 

kp  =  G^(utT^,er)  (5.12) 

where  TrG  —  0  is  required  in  view  of  (4.12).  The  validity  of  (5.12)  demands  the  same 
limitations  on  the  loading  process  as  the  validity  of  (5.10). 

Although  the  above  comments  are  sufficient  to  demonstrate  that  much  work  must 
be  done  before  a  consensus  can  be  reached,  and  that  a  formula  for  F&  based  on 
the  average  measurements  dictated  by  the  very  definition  of  cj  (e.g.,  a  commonly 
used  power  law)  is  likely  to  be  much  too  simplistic  to  encompass  all  the  phenomena 
involved,  some  concern  must  also  be  expressed  regarding  T^-dependence  of  F&. 
Indeed,  the  experiments  by  Trampczynski  et  al.  (loc.  cit.)  on  aluminum  alloys  agree 
with  a  choice  of  F ^  depending  only  on  the  second  invariant  J2  of  the  deviatoric  stress. 
In  contrast,  the  work  of  Barbee  et  al.  (loc.  cit.)  also  on  aluminum  alloys  ,  emphasizes 
a  choice  of  F &  depending  only  on  the  hydrostatic  pressure  p  =  (l/3)T,rT<^>.  Such  a 
discrepancy  can  only  be  explained  by  the  difference  in  the  experimental  procedure: 
steady  loading  in  [27]  versus  impact  in  [2].  As  a  result,  it  clearly  appears  that  the 
appropriate  form  for  F &  depends  on  the  loading  process.  At  the  least,  a  marked 
difference  must  be  made  between  quasistatic  and  quasi  instantaneous  loading.  This 
is  because  the  latter  is  accompanied  by  a  shock  wave,  the  effect  of  which  on  nucleation 
and  growth  of  microvoids  is  yet  to  be  analyzed. 

Remark  5.4:  For  all  the  purposes  enumerated  above,  one  must  bear  in  mind 
that  nucleation,  growth  and  coalescence  of  voids  within  the  grains  only  is  relevant  in 
isotropic  damage  (cf.  Remark  2.2).  This  considerably  restricts  the  relevance  of  the 
various  considerations  found  in  the  literature. 
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A  final  but  important  step  in  the  establishment  of  our  model  consists  in  checking 
its  compatibility  with  the  second  law  of  thermodynamics.  Doing  so  will  introduce 
some  limitations  in  the  choice  of  the  functions  F $  and  G ***  (assuming  that  evolution 
of  the  plastic  strain  is  governed  by  equation  (5.12)).  The  starting  point  is  the  system 
(5.7)  which  also  reads 

.  a 

Div  <r  +  p0b  =  p0u  in  £/, 

<  u  =  0  on  ro,  (5.13) 

k  <tv  =  t  on  Tj, 

upon  setting 

a  =  B(ui)[£(m)-(1-u.)-,/V1, 

6  =  6-(l/^)VrrM, 
t  =  t  +  tt(u)v. 

It  should  be  noted  that  for  x  £  U,  <r{x)  appears  to  be  the  total  stress  in  the  reference 
configuration 

a-(x)  —  T(x)  -f  Tr(u(x))I,  (5.14) 

namely  the  sum  of  the  first  Piola-Kirchhoff  and  residual  stresses  at  x.  Relation 
(5.14)  easily  follows  from  (5.5)  and  (5.1)  and  the  assumption  of  infinitesimal  elastic 
and  plastic  strains  allowing  one  to  write 

T(x)  =  (1  -  u ,(x))-2'3T^(<{>(x))  (5.15) 


upon  neglecting  higher  order  terms  in  ee  and  ep. 

Introducing  the  energy  density  xj>  =  ip(x;u,Se)  where  Se  denotes  an  arbitrary 
symmetric  3x3  tensor,  defined  by 


1 


one  finds 


0(*;u;,Se)  =  — —  B(u)S':S% 

2  po(x) 

dij> 

'dS7' 


(5.16) 


(5.17) 
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provided  that  Se  is  taken  to  be 


S’=e(u)-(  l-u>)_1/V. 


(5.18) 


In  these  notations,  the  system  (5.13)  is  readily  seen  to  be  equivalent  to  the  variational 
equality 

■ + X  ^  = X  ^ 

for  every  admissible  motion  v  (i.e.,  satisfying  v  =  0  on  r„).  Since  the  process  is 
isothermal,  the  corresponding  version  of  the  Clausius-Duhem  inequality  reads 

•  /  x  d*l>  ■  dtl> 

"  as*  -  w  a:-0' 


Equivalently,  using  (5.17)  and  (5.18),  one  finds 

cr  :  SP  -  >  0, 


with 


Sp  =  (1  -a;)-1/3ep. 


In  order  for  (5.19)  to  be  satisfied,  it  suffices  that 


<r  :  SP  >  0, 

•  dt/> 

-Wj-  >  0. 
Ou) 


(5.19) 


(5.20) 


(5.21) 

(5.22) 


From  (5.16),  an  elementary  calculation  yields 

S’)  =  {A»(TrS')2  +  %)«*„  :  S’D}  (5.23) 

where  SeD  denotes  the  deviatoric  part  of  Se  and 

K( w)  =  (1  ~u)-l/3K(u),  [,( «)  =  (1  -u,)-,/3,<M. 

with  K(u)  and  being  the  bulk  and  shear  moduli  respectively.  As  u>(x)  = 

n(0(*))  by  definition,  K(u>)  and  coincide  with  K(Q)  and  fi( 0)  of  Section  3 
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upon  substituting  fi  =  u  and  hence  are  decreasing  functions  of  u>  E  [0, 1/2].  With 
(5.23)  we  infer  that 

g<0,  0  <  u>  <  1/2, 

so  that  relation  (5.22)  will  be  satisfied  if  Co  >  0,  namely  (cf.  (5.10))  if  F^(u;,  T^)  >  0 
for  u>  E  [0,1/2].  From  (5.14)  and  (5.15),  we  may  define  the  function  F(lo,  cr)  by 

F(u,tr)  =  F^(w,T^).  (5.24) 

It  follows  that  the  condition 

F(u>,  <r)  >  0,  0  <  a;  <  1/2,  (5.25) 

is  equivalent  to  relation  (5.22).  In  words,  (5.25)  means  that  the  voids  density  u>  can 
only  increase.  Similarly,  set 

G(w,  o-,  ep)  =  G*(u>,  T^,  ep).  (5.26) 

It  is  immediate  from  (5.20)  that  condition  (5.21)  holds  if  and  only  if 
<r:(|d>ep  +  (l  -u)ep)  >  0,  0  <  a;  <  1/2. 

Recalling  (5.10)  and  (5.12)  and  from  (5.24)  and  (5.26)  above,  this  amounts  to  saying 
that 

<r:  (fF(u;,<r)ep  +  (l  -  w)G(w,  <r,  £p))  >  0,  0<u;<l/2.  (5.27) 

Thus,  with  the  choice  of  Se  as  in  (5.18)  and  u  as  internal  variables,  the  proposed 
model  is  compatible  with  the  second  law  of  thermodynamics  provided  that  inequalities 
(5.25)  and  (5.27)  hold. 

As  mentioned  earlier,  neither  relation  (5.10)  or  (5.12)  can  be  used  for  the  evolution 
of  damage  or  plastic  strain  under  nonproportional  or  cyclic  loading.  It  has  been 
known  for  long  that  invariance  of  the  plastic  strain  upon  unloading  can  conveniently 
be  handled  by  the  introduction  of  a  yield  surface  for  the  evolution  of  ep.  The  same 
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idea  can  then  be  applied  to  the  evolution  of  the  damage  u.  The  dissipative  nature  of 
the  damaging  process  suggests  the  introduction  of  a  flow  potential  P  =  P(u>,£)  where 
£  denotes  the  variable  conjugate  to  the  damage  variable,  namely  (see  also  (5.23)) 

<--£<•!«.  *->~j <5-28> 

and  Se  is  given  by  (5.18).  Evolution  of  the  damage  variable  is  then  accounted  for  by 
the  equation  (assuming  smoothness  of  the  potential  P) 

-[(8P/ae)(w,o/(dP/a «)(»,()){  if  p(u, o  =  o  and  (ap/dOMi  >  o. 

0  otherwise. 

(5.29) 

The  condition  (dP/d£)(u,  ()£  >  0  must  be  interpreted  as  a  loading  criterion.  Mono¬ 
tonicity  of  damage  thus  requires  that 

(dP/d <  0 

if  P(u>,£)  =  0  and  (dP/d£)( w,£)f  >  0.  Observing  that  both  u  and  ^  are  nonnegative 
(£  >  0  follows  from  (5.28)  and  (5.23)  and  the  growth  properties  of  the  functions  K 
and  //),  a  sufficient  condition  for  Cj  >  0  is  then  that 

(dP 

{w>0,  £>0,  P(u,Z)  =  0}  =»  (w,0  <  0, 

If  so,  the  loading  condition  reduces  to 

^>0  (5.31) 

and  the  inequality 

£u;  >  0 

ensuring  dissipation  (compliance  with  the  Clausius-Duhem  inequality)  is  automati¬ 
cally  satisfied.  The  above  formalism,  which  is  compatible  with  damage  being  unaf¬ 
fected  along  appropriate  unloading  paths,  was  suggested  by  Krajcinovic  [17]. 

It  is  rather  unlikely  that  the  potential  P  can  be  determined  on  the  basis  of  energy 
release  measurements  since  dissipation  due  to  damage  cannot  be  distinguished  from 
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dissipation  due  to  the  plastic  strain.  However,  existence  of  the  potential  P  rather 
than  its  explicit  form  is  important  for  the  theory.  Indeed,  suppose  that  P  exists  such 
that  (5.30)  is  fulfilled.  It  is  easily  checked  that  in  view  of  the  loading  criterion  (5.31), 
the  equation  (5.29)  can  be  rewritten  as 


u>t 


-[(dP/d()(u)t,(t)/(dP/duj)(uu(t)\£t  if  &  >0  and  &  =  max(max£s,£), 
0  otherwise. 


(5.32) 

where  the  subscript  “t”  relates  to  the  value  of  the  variable  at  time  t,  and  where  £ 
denotes  the  unique  real  root  of  P(0,  £)  =  0.  Uniqueness  of  £  follows  from  the  condition 
(dP/d£)( 0,£)  >  0  whenever  P(0,£)  =  0.  Existence  of  £  is  guaranteed  if  P(0,0)  <  0 
and  lim^oo  P(0,£)  >  0.  Now,  it  is  easily  seen  from  (5.16)  and  (5.17)  that  another 
expression  for  the  variable  £  introduced  in  (5.28)  is 


2Po{X) 

Therefore,  relation  (5.32)  is  equivalent  to  a  differential  equation  of  the  form 

{F{u>t,  (rt)  :  &t  if  t  >  0  and  &  =  max(max£,,£). 

«[o,fj 

0  otherwise. 

Of  course,  the  vector-valued  function  F  depends  on  the  potential  P.  But,  unlike 
P(u;,£),  the  expression  F(u,a)  :  &  directly  relates  to  void  growth  and  nucleation 
and  hence  should  be  more  easily  accessible  to  experimental  determination.  We  note 
that  similar  expressions  already  used  in  the  literature  (see  Tvergaard  [28]  and  the 
references  therein)  agree  with  the  stress  rate  &  entering  linearly.  Regarding  the 
Clausius-Duhem  inequality,  a  formula  for  F(u>,  <r)  :  &  qualifies  if  and  only  if 


F(u>t,  <Tt)  :  d’j  >  0  whenever  >  0  and  £t  =  max(max(„(), 

ae[0,t] 

where  the  number  £  also  is  to  be  determined  from  experiments. 

For  the  evolution  of  the  plastic  strain,  namely  that  of  the  variable  Sp  (cf.  (5.20)), 
the  same  procedure  applies:  the  evolution  of  Sp  may  be  governed  by  a  flow  potential 
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Q  and  the  associated  flow  rule.  Typically,  Q  =  Q(v,ir,x)  where  %  denotes  a  family 
of  suitable  hardening  parameters.  For  this,  the  reader  is  referred  to  the  standard  lit¬ 
erature  on  plasticity  theory,  except  that  the  contribution  of  damage  to  the  potential 
Q  remains  to  be  determined.  Apparently,  there  is  no  contradiction  with  basic  prin¬ 
ciples  in  introducing  two  different  potentials  P  and  Q.  This  approach  is  consistent 
with  damage  and  plastic  slip  resulting  from  distinct  phenomena  (voids  nucleation 
and  growth  versus  dislocations).  Of  course,  using  two  different  potentials  induces 
two  different  notions  for  loading  and  unloading,  but  there  is  nothing  wrong  with  that 
either:  one  merely  needs  a  “compatibility  condition”  between  P  and  Q 

ensuring  that  unloading  paths  exist  that  are  indeed  unloading  paths  for  both  P 
and  Q,  a  very  mild  restriction  in  practice. 

6  Mathematical  Aspects 

This  section  is  intended  to  show  that  the  model  for  isotropic  damage  developed  earlier 
leads  to  mathematically  tractable  boundary  value  problems.  For  simplicity,  we  shall 
limit  ourselves  to  the  quasistatic  problem  when  ti(£,  •)  is  an  equilibrium  position  for 
all  t.  This  means  that  the  forces  are  varied  sufficiently  slowly  for  the  acceleration 
term  in  (5.7)  to  be  negligible  and  the  equations  thus  take  the  form 

{Div  [jB(cu)(e(u)  —  (1  —  w)“1/3ep)J  —  V7r(u;)  +  p0b  =  0  in  U, 
u  =  0  on  ro, 

B(u)(e(u)  —  (1  —  u)~l/3ep)v  =  n(uj)v  +  t  on  Ti, 
with  evolution  laws  given  by  (5.11)  and  (5.12),  namely  (cf.  (5.24)  and  (5.26)) 

Cj  =  F(u,  <r),  u  =  0  at  t  =  0, 
kp  —  G(w,  cr.  £p),  ep  =  0  at  t  —  0. 

The  initial  conditions  are  taken  for  consistency  with  the  assumption  of  an  originally 
virgin  material  in  a  natural  configuration,  but  they  could  be  arbitrarily  chosen  for 
mathematical  purposes.  Recall  also  that  the  above  evolution  laws  are  valid  within 
the  limited  framework  of  (nearly)  proportional  and  monotone  loading  (sec  Section  5). 
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It  will  be  somewhat  more  convenient  to  use  the  variable  Sp  introduced  in  (5.20) 
instead  of  ev .  Accordingly,  we  shall  then  consider  the  problem 

f  Div  [J3(u;)(e(it)  —  Sp)]  —  V7r(u;)  +  pDb  =  0  in.  U, 


<  U  =  0  on  ro, 

(6.1) 

[  B(u>)(e(u)  —  Sp)i/  =  7r(u)i/  +  t  on  rx, 

w  =  F(u,(t), 

(6.2) 

Sv  =  ,  Sp), 

u>  =  0  at  t  =  0, 

(6.3) 

sp  =  0  at  t  =  0, 

where,  as  before 


and 


<r  =  B(w)(e(u)  -  Sp) 


(6.4) 


H(u,  <r,  S')  3  i(l  -  u )-lP(a ,  <r)S”  +  (1  -  uj)-,/3G(u>,  <r,  (1  -  u)',3Sr).  (6.5) 

In  Sections  4  and  5,  we  confined  our  attention  to  the  physically  realistic  values 
0  <  uf  <  1/2.  Here,  it  will  be  more  appropriate  to  let  u>  run  over  the  entire  real  line 
and  extend  the  functions  B,  7r,  F  and  H  by  setting 


if  u>  <  0  and 


B  M 

=  B(0)(=A(0)), 

n(u) 

=  r(0), 

F(u,a) 

=  F(0,a), 

H{u>,  o’ ,  Sp) 

=  ff(0,<7,SP), 

B{u) 

=  B(  1/2)  =  0, 

=  *-(1/2)  =  0, 

F(u,<t) 

b 

c£ 

II 

H{ui,  <r,  Sp) 

=  H(l/2,  <7 ,  Sp), 

(6.6) 


(6.7) 
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if  u  >  1  / 2. 

With  these  extensions  being  performed,  the  structure  of  the  problem  (6.1 )- (6.3) 
is  reminiscent  of  others  studied  by  Necas  and  Hlavacek  [23]  (see  also  John  [13]). 
Accordingly,  existence  and  uniqueness  results  can  be  expected  to  follow  from  the 
variational  characterization  of  u,  namely 

jv  B(ujt)e(ut)  :  e(v )  -  B(u>t)Spt  :  e(t>)  = 

/  pabt  ■  v  +  /  ir(u)t)  div  v  +  /  tt  ■  ttdA,  Vt>  £  V,  (6.8) 

v  (/  Ju  Jrx 

for  every  t  of  some  time-interval  [0,T]  (where,  again,  ujt  denotes  the  partial  mapping 
u(t,  •)  and  not  a  partial  derivative  and  similarly  with  tt(,  Sp,  etc.),  and  where  V  is 
a  suitable  space  of  admissible  displacements.  Clearly,  the  formulation  (6.8)  suggests 
the  choiceS * 7 

^  =  {v  £  (H'iU))3]  v  =  0  on  ro}  .  (6.9) 

For  consistency,  it  is  then  necessary  that  the  coefficients  of  B(ut)  belong  to  L°°(U) 
and  one  may  anticipate  that  it  will  somewhere  be  necessary  that  these  coefficients 
depend  continuously  on  u)t ■  This  condition  is  met  if  one  requires 

^  €£'([/,  n  £“(«)).  (6.10) 

But  a  difficulty  then  arises  from  the  evolution  equations  (6.2),  which,  using  (6.3), 
may  be  rewritten  as 

ut  =  f  F(u„  <r3)ds, 

Jo 

Sp  =  f  H{u>a,<r„Sp)ds- 

J  o 

Indeed,  no  matter  how  smooth  the  functions  F  and  H ,  it  follows  from  (6.4)  that 
F(u>,,cr,)  and  H(u,,  <r3,  Sp)  will  merely  be  square  integrable  since  <r3  is  no  better 
than  square  integrable  for  u,  £  V.  Of  course,  one  may  think  about  looking  for 

7Standard  notations  are  used  regarding  Sobolev  spaces. 
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u>  £  C!([/,  T],  C^{U))  instead  of  (6.10),  as  is  done  in  [13,23].  But  this  is  not  compatible 
with  the  coefficients  of  B(u>t)  being  continuous  function  of  ut  in  the  Lco(C/)-norm 
(although  these  coefficients  are  in  L°°(U)  as  soon  as  u)t  is  measurable). 

Remark  6.1:  The  above  observations  are  not  in  contradiction  with  the  results  in 
[13,23]:  a  simple  examination  reveals  that  the  existence  and  uniqueness  results  given 
there  require  much  too  stringent  assumptions  (in  many  respects)  for  our  purposes. 

To  overcome  these  difficulties,  it  is  reasonable  to  choose  the  following  modification 
of  the  original  problem:  pick  an  arbitrary  function  6  6  L°°(il3)  with  compact  support 
and  set  for  v  £  V  (e(v)  being  extended  by  0  outside  U) 

e*(v)  s  $  *  e{v)  £  C\£f)  (6.11) 

and  (compare  with  (6.4)) 

<r*  =  B(w)(e*(t»)  -  Sp).  (6.12) 

Now,  instead  of  (6.2)  and  (6.3),  prescribe  the  evolution  of  and  Sp  through 

u>  =  F(u,  tr*),  a)  =  0  at  t  =  0, 

SP  =  H(w,  <r*,  Sp),  Sp  =  0  at  t  =  0. 

Then,  the  previous  discrepancy  has  disappeared  if  F  and  H  are  continuous  and 

S- €£?«/,  n±i). 

where 

S0  =  {S  €  S,  TrS  =  0} 

and 

E  =  {5  =  (Si,)  £  (C'($)3;  S)|  =  S„}  . 

The  above  modification  of  the  original  problem  is  justified  by  the  observation  that 
e’(u)  =  e(it)  in  the  limiting  case  when  0  is  the  Dirac  delta  (so  that,  in  practice, 
0  should  be  taken  as  an  approximation  of  the  Dirac  delta).  Also,  note  the  crucial 
property  that 

v  e  v  t-*.  e*(t?)  e  s 
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is  continuous,  as  it  follows  from  standard  results  on  convolution. 


Summing  up,  we  shall  investigate  existence  and  uniqueness  of  solution  for  the 
problem 

/  B(u>t)e(ut)  :  e(v)  -  j  B{ut)Spt  :  e(v)  = 

J(J  JU 

=  J  pQbt  •  v  +  J  Tr(u>t )  div  v  +  J  tt  •  vdA,  Vv  6  V,Vt  6  [0,  T],  (6.13) 

<*=  I*  F(ua,*;)ds,  (6.14) 

Jo 

St  =  f  H(uj„  <r",  Sp)ds.  (6.15) 

Jo 

In  what  follows,  it  will  be  assumed  that  meas  ro  >  0,  F  and  H  are  locally  Lipschitz 
continuous  and 

Po  € 

6  6C'([(,co),(£6(W))3), 

t€C'([<,oo),(£e(-„))9). 

Call  uQeV  the  solution  of  (6.13)  for  t  =  0  (recall  u0  =  0  and  SPQ  =  0).  That  uQ 
exists  and  is  unique  follows  from  jB(0)  being  elliptic  and  Korn’s  inequality  (see  e.g. 
[15]).  Identify  u0  with  the  corresponding  t-independent  function  in  6],  V)  where 
S  >  0  is  arbitrary  and  let  B  denote  the  (closed)  ball  with  centre  u0  and  radius  one  in 
the  space  C*([/,6],V).  Similarly,  denote  by  B'  the  (closed)  ball  with  centre  (0,0)  and 
radius  1/4  in  the  space  8],Cl($))  x  C|([f,6],  £0)»  where  the  subscript  “o”  in  Cj 
refers  to  those  elements  vanishing  at  t  =  0. 

Pick  u  G  B.  Then,  recalling  (6.12)  and  our  assumptions  on  F  and  H ,  it  follows 
from  the  classical  theory  of  O.D.E.’s  that  the  system  (6.14)  -  (6.15)  has  a  unique 
solution 

(u>,  S')eC*([»,<I,C'(rf))  x  CftMI,  ±,). 


47 


It  is  straightforward  to  check  when  <5  >  0  is  small  enough  that  (u>,  Sp)eB'  and  that 
the  mapping 

u  £  B  *-*  (u>,  S^)  £  B',  (6.16) 

is  a  contraction  with  constant  k(S)  verifying 

lim&(6)  =  0.  (6-17) 

Conversely,  pick  (u>,Sp)  £  B'.  Since  B'  has  radius  1/4,  one  has 

B(u)t)S  :  S  >  J3(1/4)S  :  S  >  a|S|2,  (6.18) 

for  every  symmetric  3x3  tensor  S,  where  a  is  a  positive  constant  and  where 

|S|2  =  S:  S. 

It  follows  that  the  variational  problem  (6.13)  has  a  unique  solution  ut  £  V  for  every 
t  €  [0,6].  It  is  readily  checked  that  this  defines  u  as  an  element  of  Cf ([/,  6],  V)  such 
that  ut  =  uQ  for  t  —  0.  In  particular,  u  £  B  if  6  >  0  is  small  enough.  Moreover,  the 
mapping 

(w,  Sp)  £  B'  ^  u  £  B,  (6.19) 

is  Lipschitz  continuous  with  constant  C  >  0  independent  of  6.  This  is  easily  seen  from 
(6.18)  and  the  observation  that  the  mappings  B(uj)  and  7r(u>)  are  (uniformly)  Lipschitz 
continuous  tunctions  of  u>  £  R.  From  (6.17),  one  thus  finds  that  the  composition  of 
the  mappings  (6.19)  and  (6.16)  has  a  unique  fixed  point  (S  >  0  small  enough  being 
fixed).  Each  such  fixed  point  providing  a  solution  to  the  system  (6.13)-(6.15)  verifying 
the  desired  continuous  dependence  in  time,  there  follows  the  existence  and  uniqueness 
of  such  a  solution  on  [0,6]. 

Now,  let  T  >  0  denote  the  upper  bound  of  those  6’s  as  above.  Then,  a  unique 
solution  to  (6.13)-(6.15)  exists  on  [0,  T).  Suppose  that  T  <  oo  and  that 


SUP  II  IUdi=  a  <  !/2- 

6(0  ,T)  [  > 


(6.20) 
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and 


(6.21) 


Then, 


SUP  li  St  II Eo  <  °°- 

<€[0,T) 

B(ut)S  :  S  >  B(a)S  :  S  >  a|S|2,  Vt  e  [0,  T), 


(6.22) 


for  every  symmetric  3x3  tensor  S  and  some  constant  a  >  0.  Taking  v  =  ut  in 
(6.13),  it  is  immediate  from  (6.20)-(6.22)  that 


sup  jj  Ut  ||v  <  00. 
t€[0,T) 


In  turn,  with  (6.21)  and  (6.12),  the  above  yields 


sup  jj  a’  jjg  <  oo.  (6.23) 

*€[0,T) 

Combining  (6.20),  (6.21)  and  (6.23),  one  finds  that  the  functions  *)  and 

ff(cj3,  <r*,  Sp3)  are  bounded  for  s  €  [0,T).  Hence 

K  -ut2\  <  M]ti  -t2\,  ISi,  —  S£J  <  M\t\  —  *2 U  (6-24) 

for  ti,  t2  in  [0,  T),  where  M  is  a  constant  independent  of  tx  and  t2 ■  Through  Cauchy’s 
criterion,  (6.24)  shows  that 

ur  =  lim  €  C’(^), 

t—T- 

SPT  =  lim  Sp  G  E0, 

t—T~ 

exist  and  define  continuous  extensions  of  w  and  Sp  on  [0 ,T\.  Inequality  (6.22)  remains 
valid  for  t  =  T,  which  suffices  to  guarantee  the  existence  and  uniqueness  of  a  solution 
uj  £  V  of  (6.13)  for  t  =  T.  There  is  no  difficulty  in  checking  that 

uj  =  lim  ut  G  V, 

<—r- 

which  defines  a  continuous  extension  of  u  to  [0,Tj.  But  then,  it  is  possible  to  extend 
the  solution  to  the  interval  [T,  T-f  £]  for  6  >  0  small  enough  (depending  on  T )  and  this 
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extension  is  unique.  To  see  this,  it  suffices  to  repeat  the  arguments  for  existence  and 
uniqueness  on  [0,  <5]  upon  replacing  B  by  the  ball  with  centre  uj  (constant  function  of 
t)  and  radius  one  and  B'  by  the  ball  with  centre  (u>t,Sj)  and  radius  a/2  (with  a  as 
in  (6.20))  in  the  affine  subspace  of  the  space  C!([T,  T  +  S\,Cl{l£))  x  C'([T,  T  +  6],  ±,) 
of  those  elements  taking  the  value  (u>t,  Sj )  for  t  =  T.  Since  this  contradicts  the 
definition  of  T,  we  conclude  that  a  solution 

u  e  C'([/,T),  V), 
w  e  C'([/,T),C'($)), 
s"€C([',n±,), 

of  the  system  (6. 13)-(6. 15)  exists  and  is  unique  on  the  interval  [0,T)  where  T  is  the 
smallest  value  for  which  either 


SUP  II  llc.(rf)=  !/2, 

e[o,T)  1 ’ 


(6.25) 


or 

sup  ||  Spt  ||Eo=  oo.  (6.26) 

«e[o,r) 

In  view  of  (5.20),  it  is  clear  that  (6.26)  occurs  before  (6.25)  only  if  the  plastic  strain 
has  become  infinite  before  u  has  reached  the  value  1/2,  a  case  with  no  significance 
in  our  model  limited  to  small  elastic  and  plastic  strain.  If  ||  Sp  ||v0  remains  small,  so 
that  (6.25)  happens  first,  then  T  can  be  taken  as  the  failure  time  since  the  assumption 
F  >  0  guarantees  that  ut  >  0  for  every  t  £  [0,T)  and  there  must  be  points  x  £  U 
such  that  supte[0T]  ut(x)  =  1/2.  These  points  x  6  U  are  those  at  which  failure  occurs 
at  time  T  and  therefore  are  those  at  which  crack  initiation  takes  place. 


The  case  of  the  pure  traction  problem  (i.e.  T0  =  (j> )  can  be  handled  in  a  quite 
similar  way.  In  this  case,  uniqueness  is  to  be  understood  to  within  infinitesimal  rigid 
motions  and  existence  requires  the  compatibility  conditions 


/  p0bt  ■  v  —  7r(u ;()  div  v  +  /  £t  •  udA  =  0, 

Ju  Ju  Jd u 
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to  hold  for  every  t  >  0  and  every  infinitesimal  rigid  motion  v.  As  such  a  motion  is 
always  divergence-free,  these  conditions  take  the  standard  form 

/  Pobt  +  f  ttdA  =  0, 

JU  JdU 

/  p0bt  x  x  +  I  tt  x  xdA  =  0, 

JU  JdU 

for  every  t  >  0. 

For  completeness,  let  us  mention  that  in  the  case  of  the  pure  displacement  problem, 
an  alternative  proof  of  existence  and  uniqueness  of  solution  to  the  problem  (6.1  )-(6.3) 
(and  not  of  its  variational  and  modified  form  (6.13)-(6.15))  can  be  given,  at  least  if 
the  functions  F  and  H  are  smooth  enough  prior  to  be  extended  as  in  (6.6)  and  (6.7), 
such  an  extension  being  unnecessary.  Again,  the  proof  is  based  on  a  contraction 
argument,  but  solutions  are  sought  in  Holder  spaces.  More  precisely 

«€<?([/, T),  C6'“(£)), 

«eC([/,T),  r°-(£)), 

S"  eC'((i,T),  ±”'“) 

where 

T?*  =  {S  6  S1*";  TrS  =  0} 

and 

s1-”  =  (S  =  (Sij)  6  (C“"(£))3;  S>|  =  S|,}. 

Here,  a  denotes  any  real  number  such  that  0  <  a  <  1.  The  data  must  be  chosen 
accordingly,  namely 

P.  €  C‘^(£t) 

6gC'([',oo),  (C'»($)3) 

with  boundary  values  in  The  proof  no  longer  relies  on  the  variational 

approach  but  makes  essential  use  of  regularity  properties  of  the  linearized  system 
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of  elasticity,  in  particular  Theorem  6.3.7  of  Morrey  [22].  It  follows  that  while  this 
approach  can  be  applied  to  the  pure  traction  problem  as  well — assuming  t  €  C*([/,  oo), 
(C,,a(dU))*)  — it  is  not  appropriate  for  the  mixed  displacement-traction  problem. 
Recall  indeed  that  good  regularity  results  for  the  linearized  system  of  elasticity  are 
known  to  be  lacking  for  mixed  boundary  conditions.  Furthermore,  smoothness  of 
the  functions  F  and  H  seems  to  be  somewhat  atypical  and,  in  practice,  they  should 
then  be  replaced  by  smooth  approximations.  Since  this  approach  does  not  eliminate 
a  smoothing  step  of  some  kind,  the  details  of  the  corresponding  technical  proof  will 
be  omitted. 

Remark  6.2:  A  priori,  the  same  fixed  point  method  as  above  could  also  be  applied 
to  the  non-quasistatic  problem  provided  that  suitable  initial  conditions  are  prescribed 
for  u  and  du/dt.  However,  an  unexpected  difficulty  arises  from  the  fact  that  the 
mapping 

v  6  V  i->  J  7 r(u>)  div  v , 

is  obviously  not  continuous  for  the  ( L2(U ))3  -  topology.  It  follows  that  the  classical 
results  in  Lions  and  Magenes  [21]  fail  to  provide  existence  and  uniqueness  of  a  solution 
u  €  C'([/,  T],  V).  More  work  thus  needs  to  be  done  to  determine  an  appropriate  setting 
in  which  the  non-quasistatic  problem  can  be  shown  to  have  a  unique  solution  up  to 
failure. 

Remark  6.3:  If  the  evolution  laws  (6.2)  for  damage  and  plastic  strains  are  replaced 
by  suitable  flow  conditions  as  is  required  for  the  study  of  nonproportional  or  cyclic 
loading  (see  Section  5),  the  method  described  by  Korneev  and  Langer  [16]  seems  to 
be  appropriate  to  establish  existence  and  uniqueness  of  solutions  (up  to  failure)  in  a 
reasonable  approximate  problem.  Roughly  speaking,  the  approximation  consists  in 
introducing  some  inertia  in  the  yield  condition  to  eliminate  discontinuities.  But  this 
method  uses  regularity  of  the  linearized  operator  of  elasticity,  hence  is  limited  to  the 
pure  displacement  or  pure  traction  problems.  In  particular,  the  regularity  statement 


52 


in  [16,  p.  53]  about  the  mixed  problem  is  incorrect. 
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3. 


NUMERICAL  SOLUTION  OF  THE  EVOLUTION 
EQUATIONS  OF  DAMAGE  AND  RATE-DEPENDENT 

PLASTICITY 


Prologue 

This  final  chapter  is  a  reproduction  of  the  material  in  publication  #3  listed  in  Section  1.2  of 
this  report.  It  focuses  on  the  design  of  stable  schemes  for  integrating  the  evolution  equations  of 
damage  theory  and  of  viscoplasticity  which  are  robust  enough  to  model  jumps  in  the  stress  due  to 
cyclic  loading.  Extensive  numerical  experiments  are  described  and  conditions  for  the  construction 
of  effective  schemes  for  these  classes  of  problems  are  established. 

1  Introduction 

In  recent  years,  the  success  of  modeling  progressive  damage  and  rate-dependent  plasticity 
have  led  to  the  application  of  such  theories  to  an  increasing  list  of  engineering  problems.  Typically, 
such  theories  are  characterized  by  constitutive  equations  which  include  evolution  equations  for  some 
type  of  internal  variable  which  could  represent  such  features  such  as  a  loss  of  stiffness  due  to  an 
increase  in  microcrack  density,  hardness,  plastic  strain,  dislocation  density,  etc.  While  such 
phenomenological  theories  can  be  very  effective  in  modeling  history  effects,  damage,  viscoplastic 
deformation,  and  other  phenomena,  the  numerical  integration  of  the  equations  often  presents  serious 
difficulties,  particularly  when  cyclic  loading  cases  are  considered.  These  difficulties  are  related  to 
the  mathematical  stiffness  inherent  in  damage  theories  and  in  intemal-state-variable  formulations, 
with  the  result  that  many  of  the  standard  numerical  schemes,  particularly  the  explicit  schemes, 
encounter  serious  stability  or  convergence  problems. 

Several  computational  schemes  have  been  proposed  in  the  literature  for  solving 
initial-boundary-value  problems  in  viscoplasticity.  These  schemes  include  explicit  and  implicit 
methods  of  time  integration  which  are  used  in  conjunction  with  both  constant  stiffness  [1-7],  and 
tangent  stiffness  [8-11]  formulations  of  the  equilibrium  equations.  Due  to  the  stiffness  of  the 
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evolution  constitutive  equations,  many  of  these  schemes  are  only  conditionally  stable  and  may 
produce  results  which  diverge  rapidly  from  the  true  solution  when  applied  to  an  arbitrary  history  of 
loading. 

For  this  reason,  an  investigation  of  several  such  schemes  is  taken  up  in  the  present  paper  and 
two  schemes  are  identified  which  appear  to  yield  acceptable  results  when  applied  to  problems  with 
arbitrary  loading  histories.  These  techniques  include  an  Euler  forward  predictor  with  trapezoidal 
corrector  and  time  step  control,  given  in  [12],  and  Gear's  stiffy  stable  methods  with  time  steps 
selected  by  numerically  approximating  the  truncation  error,  see  Ref.  [13-14].  These  integration 
techniques  have  been  implemented  in  a  new  algorithm  which  appears  to  be  more  computationally 
efficient  than  those  previously  proposed  in  the  literature.  Our  primary  mission  here  is  to  review 
schemes  which  can  be  used  successfully  for  these  classes  of  problems  and  to  present  results  of 
applications  to  representative  problems. 

This  paper  is  divided  into  seven  sections.  Following  this  introduction,  a  brief  discussion  of 
typical  models  is  presented.  This  overview  includes  a  discussion  of  the  general  mathematical 
structure  of  these  models,  and  a  synopsis  of  several  recently  proposed  theories  of  this  type.  Section 
3  discusses  the  problem  of  mathematical  stiffness  of  the  governing  equations  and  presents  a  weak 
formulation  of  the  problem.  Section  4  and  5  are  devoted  to  the  efficiency  and  reliability  of  several 
computational  schemes  for  solving  initial-boundary-value  problems  in  intemal-state-variable 
viscoplasticity.  Here  we  are  able  to  show  that  some  methods  are  unsuitable  for  use  in  general 
purpose  finite  element  codes  while  others  appear  to  be  robust,  stable  and  efficient  for  certain 
problems.  The  final  two  sections  present  the  results  of  some  numerical  test  cases  for  two 
representative  constitutive  theories. 

2  Models  of  Viscoplasticity  and  Damage 

This  section  contains  a  brief  discussion  of  intemal-state-variable  models  for  metals  exhibiting 
time-dependent  nonelastic  deformation. 

Mathematical  Structure .  Many  of  the  intemal-state-variable  theories  and  damage 
theories  for  infinitesimal  nonelastic  deformation  have  the  same  general  structure,  typified  by  the 
following  properties: 

1.  A  strain-rate  decomposition  of  the  form 

E  =  ee  +  en 

where  e  is  the  total  strain  rate  tensor,  ee  is  the  elastic  strain  rate  tensor,  and  en  is  the  nonelastic 
strain  rate  tensor  which  includes  both  a  time  independent  inelastic  component  and  a  time  dependent 


anelastic  component. 

2.  The  nonelastic  strain  rate  is  a  function  of  the  stress,  the  damage,  and  a  set  of  internal  state 
variables 

En  =  f  (a,  d,  Zfc) 

where  a  is  the  stress  tensor,  d  is  the  damage,  and  zj.  is  a  set  of  state  variables,  and  both  d  and  zk 
may  be  tensors  and/or  scalars. 

3.  The  state  variables  vary  along  a  loading  path  according  to  certain  laws,  and  the  history 
dependence  of  the  rate  of  nonelastic  strain,  up  to  the  current  time,  is  completely  characterized  by  the 
current  values  of  the  damage  and  the  state  variables.  The  constitutive  relations  for  the  evolution  of 
the  damage  and  the  state  variables  are  of  the  general  form 

d  =  D  (a,  d,  zjj) 

Zi=gj  (a,  zk) . 

4.  Often  the  nonelastic  deformation  rate  is  deviatoric, 

tr  en  =  0. 

5.  There  need  not  exist  yield  criteria  nor  loading  or  unloading  conditions.  Hence,  nonelastic 
deformation  is  assumed  to  occur  at  all  stages  of  loading. 

These  five  properties  characterize  the  general  structure  of  most  of  the  damage  and/or 
internal-state-variable  models,  although  some  models  may  deviate  slightly  from  the  above 
descriptions.  A  list  of  some  representative  examples  follows: 


The  Bodner  and  Partom  /  Bodner  and  Stouffer’s  Theories  [15-22].  In  the  period 
1979-1983,  an  anisotropic  hardening  law  was  proposed  by  Bodner  and  Stouffer  [18,21]  in  which  a 
tensor  relationship  for  the  nonelastic  strain  rate  and  a  single  scalar  state  variable  equation  appear. 
This  model  also  uses  a  hardness  tensor  which  is  related  to  the  single  internal  state  variable  and  is 
responsible  for  the  anisotropic  material  characterization.  The  nonelastic  deformation  rate  in  an 
anisotropic  formulation  is  given  by 


£n  = 
•j 


D0Sij 


,  „  „  ,  z  ij  n  n  +1 
exp  (-0.5  ( — -  )  ( —  ) ) . 
3J2  n 


Here  Dq  is  a  scale  factor,  zjj  a  hardness  tensor,  Sy  are  the  deviatoric  stress  components,  J2  the 
second  invariant  of  the  stress  deviator,  and  n  is  a  material  constant  related  to  the  rate  sensitivity. 


57 


In  this  model,  the  single  internal  state  variable  is  the  plastic  work,  which  has  the  constitutive 


form 


Z  -  Sjjj  £jj  . 


Hart's  Theory.  [23-27].  The  equations  for  the  nonelastic  deformation  rate  are  given  by 

(a*[  II  S  -  pa  ll/p]M 

en  = -  ( S  -  4a) 

II S  -  pa  II 

where  a*,  M  and  4  are  material  constants,  S  is  the  stress  deviator,  a  is  a  tensor  internal  state 

variable,  and  II  a  II ...  Vo-  0-  .  The  evolution  equations  for  the  internal  state  variables,  o*  and  a, 
are 


a*  =  a*  Te*  /  ( ln  ( 0*  /  II  4a  I! )  )1^- 


and 


a  =  en 


[  e*/  ( ln  (a*  /  II  4a  II ) )  !A] 
II 4  a  II 


4a . 


Here  T  is  a  material  function  of  o*  and  II  4a  II,  X  is  a  material  constant,  and  £*  is  a  function 
dependent  on  the  temperature  and  a*. 

The  Gillis  &  Jones  Theory  [28].  This  theory,  proposed  for  polycrystalline  metals,  is  limited  to 
materials  having  a  linear  dependence  of  dislocation  velocity  on  the  stress  and  contains  only  one 
internal  state  variable,  the  nonelastic  strain. 

The  nonelastic  strain  rate  is 

£n  =  <}>bfp*  Vq  (p  +  0t£n)  <  0/0 Q  -  (1  +  h£n)  > 

Here  <)>  is  an  orientation  factor,  b  is  the  magnitude  of  the  Burgers  vector,  f  is  the  fraction  of  mobile 
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dislocations,  p*P  is  the  initial  dislocation  density,  p*a  is  a  dislocation  multiplication  coefficient, 
v>q  is  the  dislocation  speed  produced  by  a  stress  of  magnitude  a0,  and  h  is  a  strain  hardening 
coefficient.  In  this  expression,  <  •  >  denotes  use  of  the  step  function  <y>  =  0  if  y  <  0  and 
<y>  =  y  if  y>0. 


Robinson's  Theory  [29-30].  In  this  model,  the  multiaxial  representation  for  the  nonelastic 
deformation  rate  is  given  by 


M-l 


1  (Sy  -ay)  (Sy  -ay)  2 
En  =  -  [  - 

‘J  2  2K2 


l  (  Sjj  -OCy)  , 


where  Sjj  is  the  deviatoric  stress  tensor,  ajj  is  the  tensor-valued  internal  state  variable  (Oy=Zy) , 
and  M  and  K  are  material  constants.  The  evolutionary  equations  for  the  internal  state  variable,  given 
in  a  work  hardening-recovery  format,  are 


n 

2(lH  £y 

ay  ay  P/2 
2K2 


aklakl 

R - 

2K2 


n-JLl 

2 

“ij  * 


Here  |i,  H,  n,  p  and  R  are  also  material  constants. 


Miller's  Theory  [31-32].  The  development  of  this  model  stems  from  a  desire  to  accurately  model 
steady-state  creep  deformation.  Using  this  as  a  basis,  the  following  one-dimensional  generalized 
form  for  the  nonelastic  deformation  rate  is  obtained: 


where 


en  =  B8'  [  sinh  ( 


lo-RI 


1.5  n 


)] 


sgn  (o-R)  , 
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8'  =  e*p  [(06i^ :)  ‘  (ln<_T -  )  +  1)]  ;  TS0'6Tm 

=  e*P  [  ;  T  >  0.6Tm 

In  these  expressions  B,  and  n  are  material  constants,  R  and  D  are  internal  state  variables 
representing  the  back  stress  and  drag  stress,  respectively,  a  is  the  applied  stress,  Q  is  an  activation 
energy,  K  is  Boltzman's  constant,  T  is  the  temperature,  Tm  is  the  melting  temperature,  and  sgn 
represents  the  signum  function.  The  evolution  equations  for  the  internal  state  variables  are 

R  =  H!  en  -  Hj  B0'  (sinh  (A!  IRI )  )n  sgn  (R) 

D  =  H2  l£nl  (C2  +  IRI  -  A2  /  A! )  D3  -  H2  C2B0'  (sinh  (A2D3))n  . 

Here  Hj,  H2,  C2,  Aj,  and  A2  are  also  material  constants. 


The  Krieg,  Swearengen  and  Rhodes  Theory  [33].  In  1978,  Krieg  et  al.  proposed  a  unified 
model  for  creep  and  plasticity  in  metals,  using  two  internal  state  variables  to  reflect  the  current 
microstructural  state.  The  model  is  applicable  to  metals  under  isothermal  conditions  in  the 
temperature  range  of  0.3  Tm  to  0.7  Tm,  where  Tm  is  the  homologuous  temperature. 

The  constitutive  equations  for  the  non-elastic  deformation  rate  are  based  on  a  general  function 
form  for  dislocation  glide  given  by 


En 


IIS  -  a|l  -,m 
R  -* 


S  -  a 
IIS  -  all  ' 


In  this  expression,  S  is  the  deviatoric  applied  stress,  a  is  the  back  stress  tensor,  R  is  the  drag 
stress,  £g  and  m  are  temperature-dependent  material  constants,  and  llo||  =  Vdy  cty.  The  evolution 
equations  for  the  two  internal  state  variables  are  postulated  as 


raa 

a  =  AorE11  - -  and  R  =  Ad  II  en  II  -  ra 

Hall  R  a 
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where  Aa,  Ar,  ra  and  rR  are  hardening  and  recovery  functions,  respectively;.  Specific  forms  for 
these  hardening  and  recovery  functions  are  given  in  [33]. 


Cernecky  and  Krempl's  Theory  [34-38].  In  1980,  Krempl  and  associates  proposed  a  coupled 
infinitesimal  theory  of  thermo-visco-plasticity.  This  constitutive  model  is  based  on  a  nonlinear, 
multiaxial  generalization  of  the  standard  linear  solid,  which  consists  of  a  spring  in  parallel  with  a 
Maxwell  element  This  generalization  leads  to  a  proposed  constitutive  form  of 

^pq  "*■  Kpqmn  (O,  £,  T)  (7mn  =  Gpq  (£,  T)  +  Mpqmn  (O',  £,  T)  (2.1) 

where  Kpqmn ,  G^,  and  Mpqmn  are  material  functions  determined  from  experimental  data,  and 

Kijmn  Mmnkl  =  Djjid  (2.2) 


DjjRi  being  the  tensor  of  linear  elastic  moduli.  Solving  equation  (2.1)  for  amn  and  using  equation 
(2.2)  gives  the  following  relation  for  the  nonelastic  strain  rate 

Ey  =  Mjjmn  (o,  £,  T)  [<Jmn  -  Gmn  (e,  T)  ]  (2.3) 

Also  by  using  the  relation 

-1  n 

emn  =  ^ijmn  °kl  +  emn 

for  the  total  strain  rate  in  (2.3),  a  final  form  for  the  nonelastic  strain  rate  is  obtained  which  is  a 
function  only  of  the  stress  and  internal  variable. 


The  Cescotto  and  Leckie  Theory  [39].  The  constitutive  forms  for  the  nonelastic  strain  rate 
and  internal  variables  are 


a-a2 

€n  =  f  ( I - 1 )  sgn  (0-0.2) 

a3 


61 


en-raa2 


and 

a3  =  hk  -  rk 

where  f  is  the  nonelastic  strain  rate  function,  ha  and  hk  are  hardening  functions,  and  ha  and  rk 
are  recovery  functions.  Here  sgn  represents  the  signum  function  which  takes  on  the  values  of  1, 0, 
and  -1,  depending  on  whether  the  argument  is  positive,  zero  or  negative.  No  particular  form  of  the 
functions  f ,  ha ,  hk ,  ra  and  rk  is  assumed  a  priori,  and  a  set  of  experiments  is  required  to 
define  those  functions. 


3  Stiffness  and  Approximation 

A  wide  range  of  techniques  has  been  proposed  in  the  literature  to  numerically  integrate  the 
evolution  equations  described  above.  These  methods  include  constant  stiffness  and  tangent  stiffness 
formulations  of  the  equilibrium  equations  used  in  conjunction  with  either  explicit  or  implicit  time 
integration  techniques.  Testing  of  these  techniques  is  often  done  in  the  context  of  a  monotonic 
loading  situation  and  for  restricted  forms  of  the  constitutive  equations,  even  though  arbitrary 
loading  histories,  especially  cyclic  loading  situations  often  present  very  severe  numerical 
difficulties. 

In  this  section,  we  examine  several  computational  methods  and  integration  techniques.  We 
focus  on  computational  issues,  stability  and  overall  performance  of  representative  algorithms.  For 
definiteness,  we  confine  our  attention  to  internal  state  variable  theories  appropriate  for  quasi-static, 
infinitesimal,  isothermal  deformations. 

The  general  form  of  the  constitutive  relationships  given  above  is: 


en  =  f(o,  4) 

Zi  =  gi  (<*>  Zk) 

a=E(e-e")  =  E  ( e  -  f  (a,  Zk) ) . 
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where  here  can  also  represent  a  damage  measure.  These  differential  equations  are  to  be  integrated 
over  a  time  interval  (0,T)  subject  to  the  constraint  imposed  by  the  equilibrium  equations  on  the  total 

strain  rate.  Assuming  that  the  total  strain  rates  are  prescribed  over  time,  it  is  obvious  that  Aen  is 
available  upon  inverting 

Ao  =  E  ( Ae  -  A£n ) . 

Thus  it  is  desirable  to  integrate  a  and  Zj  subject  to  the  equilibrium  constraint. 

If  e  is  prescribed,  then  the  above  problem  reduces  to  a  system  of  the  form 


o  =  k  (a,  Zk) 

=  Si  (°»  Zk) 


which  can  be  rewritten  simply  as  the  dynamical  system, 

y  =  F  (y) . 

Such  a  system  of  differential  equations  given  by  y  =  F  (y)  is  said  to  be  stiff  (see  Lambert 

[40]  if 

(1)  <  0  for  i  =  1,2, ...  m 

(2)  max  I  Re  Xj  I  »  min  I  Re  I 

i  i 

where  are  the  eigenvalues  of  the  Jacobian  matrix  3F  /  dy  and  m  is  the  number  of  equations 
in  the  system. 

Thus,  a  system  of  stiff  differential  equations  is  one  in  which  the  components  of  the  solution 
may  be  changing  or  decaying  at  greatly  different  rates  over  a  time  interval.  Then  the  evolution  of  the 
rapidly  decaying  component  of  the  solution  may  require  very  small  time  steps  to  be  used  in  a 
numerical  scheme  while  the  component  associated  with  the  largest  eigenvalue  may  need  to  be 
integrated  over  a  relatively  large  period  of  time. 

In  the  classical  numerical  solutions  to  ordinary  differential  equations,  the  problem  of 
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numerical  stability  concerns  the  growth  of  truncation  and  round-off  errors  from  one  time  step  to 
another  in  a  given  numerical  integration  scheme.  A  numerical  integration  scheme  is  said  to  be 
absolutely  stable  for  a  given  time  step  At  and  differential  equation,  if  the  change  in  the  solution  due 
to  a  perturbation  5  in  one  mesh  value  yn  is  no  larger  than  5  in  all  subsequent  values  ym  for  n  < 
m  (cf.  Gear  [41]). 

For  the  standard  test  problem,  x  =  Xx,  X  being  constant,  one  can  define  the  region  of 
absolute  stability  of  a  given  algorithm  as  that  set  of  values  At  and  X  for  which  a  perturbation  in  a 
single  value  xn  will  produce  subsequent  values  which  do  not  increase  from  time  step  to  time  step. 
Thus,  for  the  forward  Euler  integration  of  the  test  equation  x  =  Xx  given  by 

xn+l  =  xn+At  xn  =  (1  +  X  At)  xn  . 

there  is  no  region  of  absolute  stability  for  X  >  0.  For  X<  0  it  is  necessary  that  1 1+XAt  Id  so  that 
XAt  lies  in  the  unit  circle  centered  at  -1  in  the  XAt  complex  plane.  For  these  values  of  XAt  inside 
the  unit  circle,  the  integration  may  be  performed  without  errors  growing  from  one  time  step  to 
another.  Similarly,  regions  of  absolute  stability  may  be  determined  for  any  numerical  method. 

Sample  Time  Step  Calculations.  For  the  one-dimensional  form  of  the  constitutive  equations 
of  Bodner  et  al.,  eigenvalue  calculations  have  been  performed  for  various  materials  and  total  strain 
rates.  These  eigenvalues  were  then  used  to  estimate  stable  time  steps  for  the  forward  Euler  type 
integration.  The  results  are  presented  in  Figs.  1  and  2  with  the  stable  time  steps  set  to  0.65  seconds 
in  the  initital  "elastic  region"  where,  otherwise,  very  large  values  would  have  been  obtained.  From 
this  data  large  variations  in  eigenvalues  can  be  seen  for  different  regions  of  strain,  for  various 
materials  and  strain  rates.  It  is  notable  that  for  the  strain  rate  of  2.0  E-4  /  sec.  to  a  strain  of  2 
percent,  say,  may  require  between  12,000  and  20,000  time  steps  for  the  copper  or  aluminum 
specimens. 

The  shapes  of  these  time  step  curves  may  also  be  used  to  explain  the  oscillations  in  some 
numerically  obtained  stress-strain  curves  such  as  that  in  Fig.  3.  In  this  figure,  time  steps  were 
initially  selected  slightly  larger  than  the  acceptable  stable  time  step.  Therefore,  on  exiting  the  elastic 
region,  errors  were  introduced  which  were  oscillatory  in  nature  but  not  catastrophic.  As  the 
integration  continued,  the  stable  time  step  increased  and  the  oscillations  decayed  as  the  time  step 
entered  the  stable  region. 

Weak  Formulation.  Subsequent  calculations  are  performed  on  systems  resulting  from  a  finite 
element  approximation  of  weak  forms  of  the  momer  m  equations. 
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where 


e 

(Eijkl  eJd)j  =  -  bi 


in  £2 


e  n 

eki  =  U(k,j)  -  £k] 


Ekl  ~  ^kl  (®»  Zfc  ) 


(3.1) 


(3.2) 


Zj  =  gi  (O,  Zk) . 


Here 


a  =  JCadt+o,  o  =  E  ee  (3.3) 

0 

Zi  =  JtZidt  +  Zi  .  (3.4) 

0 

A  weak  form  of  (3.1)  is  obtained  in  the  usual  way:  multiply  (3.1)  by  a  suitably  smooth  test 
function  o,  =  u^x),  integrate  over  Q,  and  use  the  Green-Gauss  divergence  theorem  to  integrate 
by  parts  the  stress-power  terms.  Let  V  denote  the  space  of  test  functions 

V  =  {  Uj  G  Wm-P  (Q)  1 =  0  a.e.  on  9Qj,  1  <  i  <  N  } 

where  Wm’P(Q)  is  the  Sobolev  space  of  order  (m,p),  with  m>0,  meIR,  l<p<°o  and  where 
specific  values  of  m  and  p  depend  upon  the  particular  forms  of  the  constitutive  equations 
governing  the  material  under  consideration,  (for  the  cases  considered  here  m  =  1 ,  p  =  2).  The 
weak  form  of  the  boundary-initial-value  problem  (3.1)  -  (3.4)  is  then: 


Find  a  displacement  rate  field  t  1-4  u  (x,t)  €  V+  {u }  such  that  for  every  t  e  [0,T], 

J  Ejjjd  ukJ  vi  j  ^  =  J  ^ijlcl  £kl  (u)  j 

ft  ft 

+  JbiuidQ+J  TjUjds  V  \),  e  V  (3.5) 
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ft  an2 

with 

n 

£kl  =  fkl  (a*  Zk) 

and  a  and  Zj  are  given  by  (3.3)  and  (3.4).  Here  dO.  and  ds  are  volume  and  surface  measures, 
u  is  any  function  defined  on  Q  such  that  its  trace  on  the  boundary  segment  is  u  ( where 

ui  I  0q^=  uj),  and  a  and  Zj  are  understood  to  depend  upon  ujj  through  (3.3)  and  (3.4),  and  to 

satisfy  initial  conditions  CTjj  (x,  0)  =  j  (x),  Zj  (x,  0)  =  Zj  (x).  It  is  easily  verified  that  any 
sufficiently  smooth  solution  of  (3.5)  will  also  satisfy  the  governing  equations.  Conversely,  any 
solution  of  the  governing  equations  and  boundary  conditions  will  also  satisfy  (3.5). 

Using  standard  notations,  a  finite  element  approximation  of  (3.5)  leads  to  the  discrete  system 
of  evolution  equations.  If  the  discrete  velocity  components  are  of  the  form 

i  N  i 

Uh(X,t)  =  X  Uj(t)<j>j(x) 

j=l  (3.6) 

where  N  is  the  number  of  nodes,  i  indicates  the  vector  component,  <j)j  is  a  basis  function,  Dj  is 
the  value  of  the  test  function  at  node  xj ,  and  uj(t)  is  the  value  of  u^  at  node  xj  at  time  t, 
then  we  wish  to  find  the  vector  of  nodel  displacements  u  such  that 

f  BT  E  B  u  dQ  =  f  BT  E  en  dQ  +  f  d>T  b  dC2  +  f  <J>T  T  ds 

ft  ft  Jft  dn2 

(en  =  f(o(u),Zk)  ,  Zk  =  g(a,Zk)  .  (3.7) 

4  Algorithms  for  Integrating  Stiff  Evolution  Equations  of  Viscoplasticity 

We  now  outline  three  popular  algorithms  for  rate-dependent  plasticity  found  in  the  literature 
and  propose  an  alternative  scheme  that  performs  very  well  for  cyclic  loading  cases. 

The  general  strategy  in  these  algorithms  is  as  follows:  with  the  initial  distribution  of  the  stress 
and  internal  variables  specified,  use  the  equilibrium  equations  to  supply  the  spatial  variation  of  the 
constraint  (the  momentum  equations).  Then  integrate  the  constitutive  equations  forward  in  time, 
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evaluating  principal  variables  at  the  integration  points.  With  the  updated  values  of  the  stress  and 
internal  variables  at  the  new  time,  the  constraint  condition  is  again  imposed.  This  sequence  of 
determining  the  constraint,  then  advancing  the  constitutive  equations  in  time  is  continued  until  the 
desired  history  of  the  initial-boundary-value  problem  has  been  traced. 

The  Initial  Strain  Rate  Method.  This  method  was  originally  proposed  in  1972  in  [7].  Starting 
with  governing  differential  equations  in  the  rate  form,  a  finite  element  approximation  of  the 
equilibrium  equations  is  constructed  (as  in  (3.7)),  giving 

J  BT  E  B  u  dQ  =  [  Bt  E  en  dQ  +  F  (4.1) 

ft  ft 

where  F  is  the  vector  of  force  rates  derived  from  surface  traction  rates  and  body  force  rates.  The 
algorithm  is  then: 

1 .  Initialize  ct  ;  Zj ,  set  tn  =  0 . 

2.  Calculate  en  =  f  (ct,  Zk)  at  t  =  tn  (thus  determining  the  right-hand-side  of  (4.1) ). 

3 .  Solve  the  equilibrium  condition  for  %  ( t„) . 

4.  Calculate  e(tn)  =  B  uh  ( tn) . 

5.  Calculate  o(tn)  =  E  (e(tn)  -  en(tn) ) . 

6.  Calculate  Zj  (tn)  =  gj  (o(  tn) ,  Zk  (tn) ) . 

7.  Integrate  ct,  Zj  forward  over  some  appropriate  At  to  CT(tn+1),  Zj(tn+j). 

8.  If  tn  =  At  <  T  to  (2) . 

9.  Stop. 

Step  (7)  characterizes  an  explicit  scheme,  and  this  step  can  be  displaced  by  a  subroutine  for  implicit 
method,  if  needed.  If  predictor-corrector  type  integration  schemes  are  used,  a  slight  modification  is 
necessary  beginning  with  (7). 

0  0 

7.  Estimate  ct  (tn+i),  Zj  (tn+1)  using  a  predictor. 

8.  Solve  (4.1)  for  uj,  (tn+1)  using  latest  entries  for  c(tn+i),  Zj  (tn+i) . 

9.  Calculate  e'On+i)  and  a'(tn+i)- 

c 

10.  Use  a  corrector  to  update  o^tn+i),  Z,  (tn+j) . 

11.  Check  an  appropriate  tolerance  index  for  convergence.  If  not  achieved,  go  to  (8)  with  a 
new  estimate  of  CT(tn+1)  and  Zj  (tn+1) . 


12.  Set  o(tn+i)  —  <Atn+i)  snd  ^iOn+i)  On+l)  • 

13.  If  t  + At<T  to  2 . 

14.  Stop. 

The  Initial  Strain  Method.  This  method  differs  from  that  above  in  that  the  equilibrium 
equations  are  written  in  incremental  form  rather  than  rate  form.  This  results  in  a  finite  element 
approximation  of  the  equilibrium  equations  of  the  form 

/  BT  E  B  Au  d£2  =  J  BT  E  Aen  dfl  +  AF  (4.2) 

ci  ci 

where  A  represents  an  increment  of  the  proposed  quantity. 

Thus  the  computational  scheme  becomes  (for  the  case  of  explicit  integration  in  steps  (2)  and 
(7)): 

1.  Initialize  cr,  Zj ,  set  tn  =  0 . 

2.  Calculate  Aen  over  At . 

3.  Solve  (4.2)  for  A  iq, . 

4.  Calculate  Ae  over  At  from  Ae  =  BA  uj, . 

5.  Calculate  Ao=  E  (Ae  -  Aen ) . 

6.  Calculate  Zj  =  gi(o(tn),  Zk(tn)). 

7.  Integrate  Z,  forward  over  At  to  Z4  (tn+i) . 

8.  Set  o(tn+1)  =  a(tn)  +  A  a. 

9.  If  tn  +  At  <  T  to  (2) . 

10.  Stop. 

A  similar  revision  to  that  discussed  earlier  for  implementation  of  predictor-corrector  methods  can  be 
made. 

A  Forward  Gradient  Scheme.  This  method  is  considerably  different  from  the  methods 
introduced  above.  It  involves  a  tangent  stiffness  matrix  calculation  updated  from  step  to  step,  and  a 
particular  type  of  numerical  integration  specified  for  the  time  integration. 

Beginning  with  the  incremental  form  of  the  equilibrium  equations,  we  have 

J  BTEBAudQ=J  BTEAendn  +  AF 
n  a 


or 
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J  BTAodn  =  AF. 
a 

It  is  assumed  that  Aen  and  AZj  are  given  by 

Aen  =  At  {  (1  -  0)  E"(tn)  +  0£"(tn+1)} 

AZj  =  At  {  (1-0)  Zi(tn)  +  eZi  (tn+1)}. 

Expanding  en  (tn+i)  and  Zj  (tn+1)  in  a  Taylor  series  about  tn  results  in 

Ae"  =  At  {  en  (t„)  +  0  An  Aa  +  6  Bn  AZj } 

AZj  =  At  {  Zj  (tn)  +  0  Q,  Aa+  0  D„  AZj} 

where 

3  £n  3en  3Zj  3Zj 

An  =  -^Y :  L  *  =  -jyjr  *  »  Qi  = '-TT  ]  »  =  I  • 

oa  t„  dZj  tn  do  t„  dZk  tn 


(4.3) 

(4.4) 


Solving  for  AZj  gives 

AZj  =  [  I  -  0At  Dn]-1  At  [  Zj  (O  +  0Cn  Aa]  .  (4.5) 

Substituting  this  into  (4.3)  for  AZj  and  neglecting  terms  of  order  At2  gives 

Ae"  =  £n  (tn)At  +  0AtAn  Ao  . 

Thus, 

Ao  =  E  [Ae  -  Ae"]  =  [  I  +  0At  EAj1  E[BAu  -  At£n(tn )] 

=  D  [  BAu  -  At  £n(tn )]  .  (4.6) 

Finally,  substituting  this  into  the  equilibrium  equation  gives 
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The  computational  algorithm  then  becomes 


1. Initialize  a;  Z*,set  tn  =  0. 

2.  Calculate  en  (tn  )  =  f  (a  (tn ),  Zj  (tn )) . 

3.  Solve  (4.7)  for  A  uj, . 

4.  Calculate  Ae  =  BA  uh . 

5.  Calculate  A  a  from  (4.6) . 

6.  Calculate  AZj  from  (4.5) . 

7.  Update  o(tn+i)  =  a  (^ )  +  Aa 

Zj  (Wl)  =  Zj  (tn)  +  A  Zj . 

8.  If  t  +  At  <  T  set  tn  =  t  +  At  go  to  (2) . 

10.  Stop. 

New  Algorithm.  A  new  computational  method  suggests  itself,  which  is  similar  to  the  initial 
strain  algorithm  in  that  an  incremental  form  of  the  equilibrium  constraint  condition  is  imposed  as 

f  BT E  B  Au  dQ  =  f  BTEAenda  +AF.  (4.8) 

QQ 

We  then  proceed  as  follows: 

1 .  Initialize  o ,  Zj ,  set  tn  =  0 ,  and  select  AT . 

2.  Estimate  Ae„  over  AT  (using  predictor  method) . 

3.  Solve  the  equilibrium  constraint  for  Auj,. 

4.  Calculate  Ae  =  BAu^. 

5.  Assume  e  =  Ae/AT  is  constant  over  AT. 

6.  Subincrement  and  integrate  a,  Zj  accurately  over  AT,  neglecting  the  previous  estimate 
of  Ae"  . 

7.  Calculate  a  new  guess  of  Ae"  =  Ae  -  E'1  Aa. 

8.  Check  for  convergence  of  Ae0  with  Ae0;  if  no  convergence  occurs  set  A£n  =  A£n  and 
go  to  (3),  otherwise  set  Ag  =  E(A£-Ae°)  . 

9.  If  tn  +  AT  =  T ,  Stop 

10.  Select  a  new  AT  and  go  to  (2). 

This  method  possesses  the  following  desirable  characteristics: 


1 .  A  constant  elastic  stiffness  matrix  is  used  throughout  the  solution  process. 

2.  Different  time  steps  at  different  integration  points  allowed  and  the  overall  AT  is  not 
restricted  to  be  the  smallest  time  step. 

3.  Integration  using  subincrement  requires  the  coordinates  and  geometry  of  each  elment  to  be 
loaded  much  less  frequently  in  the  finite  element  simulation. 

4.  By  selecting  a  larger  overall  time  step  AT,  fewer  total  time  steps  may  be  taken  and  thus 
fewer  right  hand  sides  need  be  considered. 

5  Numerical  Integration 

This  section  presents  techniques  for  integration  of  the  constitutive  equations  which  constitutes  a 
critical  step  in  the  algorithms  discussed  in  Section  4.  To  test  the  applicability  of  various  integration 
techniques,  a  group  of  one-dimensional  test  problems  has  been  selected  which  covers  a  variety  of 
loading  histories.  These  problems  have  been  taken  from  examples  in  the  literature  to  ensure  that  the 
constitutive  models  are  used  in  the  correct  context.  The  constitutive  models  being  integrated  are 
those  of  Bodner  and  associates  and  Hart  which  are  representative  of  the  general  exponential  or 
power  law  form  often  assumed  for  the  nonelastic  strain  rates. 

The  first  technique  to  be  considered  is  the  simple  forward  Euler  method.  This  method  is  the 
simplest  and  easiest  to  use  but  suffers  from  being  a  first-order  method  and  only  conditionally  stable. 
To  use  this  method,  it  is  necessary  to  select  time  steps  so  that  the  CFL  conditions  on  stability  are  not 
violated  and  an  acceptably  small  truncation  error  is  introduced.  This  selection  of  a  suitable  time  step 
for  arbitrary  loading  histories,  constitutive  equations,  and  material  types  prohibits  the  method  from 
being  useful  without  some  kind  of  automatic  time  step  selection. 

A  first  possibility  of  time  step  selection  is  that  of  calculating  the  maximum  eigenvalue  at  each 
point  in  time.  This  calculation  is,  however,  too  time  consuming  and  is  also  generally  unacceptable 
since  the  maximum  eignevalue  may  be  changing  over  the  time  step. 

A  second  possibility,  proposed  in  [2,3],  consists  of  selecting  a  time  step  so  that  the  increment 
of  nonelastic  strain  over  the  step  is  some  fraction  of  the  total  strain.  Thus,  a  proposed  step  is  given 
by  At  =  xe/en,  where  x  is  a  constant,  generally  around  0.1,  and  the  bars  represent  equivalent 
values  of  the  indicated  quantities.  Unfortunately,  this  method  of  time  step  selection  is  also  generally 
unacceptable  because,  for  stress  states  in  the  "elastic  region",  en  =  0  and  a  very  large  time  step  is 
them  indicated.  Results  obtained  using  this  type  of  step  selection,  with  a  maximum  allowable  step 
size  imposed,  are  similar  to  those  shown  in  Fig.  4. 

A  third  method  of  time  step  selection  is  presented  in  [12]  for  a  single  differential  equation 
y  =  F(y).  This  time  step  control  is  based  on  a  comparison  of  a  suitably  defined  error 
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I  At*  (y(tk)  -  y(Vi) )  I 
I  y(tk) 1 


with  prescribed  error  limits  emax  and  emjn .  The  time  step  at  the  k-th  step,  Atk ,  is  then  defined 
on  the  basis  of  its  estimate  Atk ,  according  to 

emax<e  :  replace  At^  by  At^  and  recompute  e 

e  -emax  :  setAtk  =  At^  and  compute  yC^+j) 

where  the  initial  time  step  is  prescribed.  The  next  step  size  is  then  estimated  by 

emin <  e  —  emax  •  se*  ^k+l  = 
e  <  emin  •  se*  ^k+1  =  2Atk. 

An  extension  of  this  technique  to  the  vector  case  y  =  F(y)  is  accomplished  by  introducing  the 
infinity  error  norm  and  thus  selecting  the  maximum  value  of  e  from  all  vector  components  to 
determine  an  acceptable  step  size. 

This  a-priori  method  of  time  step  selection  suffers  from  deficiencies  similar  to  those  of  the 
second  time  step  method,  as  can  be  seen  from  Fig.  5.  Satisfactory  results  were  obtained,  however, 
for  the  constraint  strain  rate  simulations  (see  Figs.  7  - 1 1)  the  strain  rate  change  tests  (Figs.  12-13) 
and  the  stress  relaxation  tests  (see  Fig.  15).  Results  obtained  by  using  this  method  on  the  loading, 
unloading  and  reloading  test,  and  also  in  cyclic  testing  situations,  as  is  demonstrated  in  Fig.  5,  were 
less  than  satisfactory.  Note  that  the  success  of  using  this  method  depends  strongly  on  the  size  of 
the  initial  step  size  prescribed. 

For  initial  steps  that  are  "too  large",  stability  difficulties  or  spurious  loading  paths  may  be 
encountered  in  any  of  the  test  problems.  Thus  it  appears  that  the  forward  Euler  method  with  or 
without  time  step  control  is  probably  not  a  good  choice. 

Another  group  of  integration  techniques  suggested  in  the  literature  are  the  higher-order  explicit 
schemes.  These  techniques  which  include  second-  and  fourth-order  Runge-Kutta,  and 
Adams-Bashford  methods,  and  are  suggested  by  some  authors  with  the  hope  that  higher  order 
methods  will  allow  larger  time  steps  because  of  thier  inherently  smaller  truncation  errors. 
Unfortuantely,  for  such  methods  the  issue  of  conditional  stability  is  of  central  importance  and 
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dominates  the  issue  of  time-step  selection.  These  methods  have  regions  of  absolute  stability  which 
are  approximately  equivalent  to  the  stability  region  of  the  forward  Euler  method,  with  the  result  that 
often  no  appreciable  difference  can  be  seen  between  the  results  of  these  higher  order  methods  and 
the  forward  Eule  method. 

A  third  group  of  integration  techniques  are  the  implicit  schemes  which  employ  direct  or  Jacobi 
iteration  in  the  solution  of  the  resulting  nonlinear  equations.  These  methods  may  be  expressed,  for 
a  single  equation  y  =  F(y),  in  a  general  form  as 

yn+i  =  PAtF(yn+i)  +  q  (5.  i> 

where  n  +  1  indicates  evaluation  at  time  tn+j ,  At  is  the  time  step,  P  is  a  constant  dependent  on 
the  numerical  scheme,  and  q  is  a  known  function  of  eviously  calculated  values  of  y  and  F. 
Using  direct  or  Jacobi  iteration  implies  that 

yn+l  =  PAt  F(yn+1)  +  q  =  $  (yn+i)  (5.2) 

and  (s)  indicating  the  iteration  number  and  yn+ j  is  provided  by  some  predictor  calculation.  This 
sequence  of  approximations  given  by  (5.2)  will  converge  to  the  solution  of  (5.1)  whenever 
4>  (yn+i)  satisfies  a  Lipschitz  condition 

*  * 

•Wyn+l)  -  <>(yn+l)l<  Mlyn+1  -  yn+I  I 
* 

for  all  yn+j  and  yn+j ,  where  the  Lipschitz  constant  M  satisfies  0  <  M  <  1.  Then  there  exists  a 

unique  solution  y  of  (5.1),  and  the  sequence  of  approximations  defined  by  (5.2)  is  such  that 

lim  y(s)  =  y  (see  [40]) .  A  similar  result  also  applies  for  systems  of  equations  with  the  absolute 

values  replaced  by  norms  of  corresponding  vectors. 

Integration  techniques  suggested  in  the  literature  which  fall  into  this  category  include  forward 
Euler  predictor  with  trapezoidal  corrector,  forward  Euler  predictor  with  backward  Euler  corrector, 
and  mid-step  integration  where  rates  at  the  midpoint  of  the  time  step  are  used  in  the  integration. 
These  methods  are  superior  to  those  discussed  in  the  proceeding  paragraphs,  in  that  they  have  large 
regions  of  absolute  stability,  which  include  the  negative  complex  half  plane  and  essentially  eliminate 
the  stiffness  difficulty.  They  are,  however,  limited  by  the  radius  of  convergence  of  the  direct 
iteration  technique.  For  example,  if  the  partial  derivatives  of  the  Jacobian  matrix  J  =  dF/dy  are 
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continuous  and  bounded  in  an  appropriate  region,  then  the  Lipschitz  constant  of  F  may  be 

taken  to  L  =  II J II.  Now  for  any  matrix  A,  we  have  II A  II  £  p  (A)  where  p  (A)  =  max  I  Xi  I 

i 

and  Xi  are  the  eigenvalues  of  J.  This  implies  that  the  Lipschitz  constant  L,  is  greater  than  or  equal 
to  the  magnitude  of  the  maximum  eigenvalue.  Considering  the  general  form  given  in  (5.1),  we  can 
choose  the  Lipschitz  constant  M  to  be  LAt  l(3|  which  implies  that  (5.2)  converges  to  the  solution  of 
(5.1)  if  At  <  1/  (Llpi)  or  approximately  At  <  l/(ipiX.max)-  This  result  suggests  a  time  step 
restriction,  again  related  to  the  maximum  eigenvalue,  which  may  be  as  restrictive  as  the  stability 
requirements  of  the  previous  methods. 

If  these  predictor  corrector  methods  with  direct  iteration  are  to  be  used,  some  type  of  time  step 
control  is  again  necessary.  One  possible  time  step  selection  technique  for  use  with  the  forward 
Euler  predictor,  trapezoidal  corrector  is  given  in  [12]  for  a  single  differential  equation  y  =  F(y). 
This  automaic  control  is  implicit  in  nature  in  that  an  initial  time  step  size  is  estimated  from  the 
previous  time  step  and  then  adjusted  after  the  corrector  calculations  have  been  performed.  This  time 
step  control  is  basically  the  same  as  the  third  technique  described  above  for  use  with  the  forward 
Euler  method,  except  that  for  this  case  the  error  is  defined  by 

IAtk  (Ffc+1  -  Fk)l 
e  =  - 

2iyk+i* 


c  p  p 

where  yk+i  is  the  corrector  value  of  y^+j  and  Ffc+j  =  F(yjc+j)  . 

This  combination  of  predictor-corrector  with  time  step  control  performed  satisfactorily  on  all 
the  test  problems,  but  was  not  efficient  in  test  situations  when  the  constitutive  equations  were  fairly 
stiff.  Results  for  the  constant  strain  rate  simulations,  the  strain  rate  change  test,  the  stress  relaxation 
tests,  creep  test,  and  the  stress  change  test  are  the  same  as  those  obtained  by  the  forward  Euler 
method  with  time  step  control,  with  plotting  accuracy.  Use  of  this  combination  on  the  loading, 
unloading,  reloading  test  is  shown  in  Fig.  14,  on  the  cyclic  tension  compression  tests  is  shown  in 
Figs.  19  -  22,  on  the  cyclic  relaxation  test  in  Fig.  23,  and  on  the  cyclic  creep  test  in  Fig.  24. 

We  also  mention  a  very  popular  group  of  integration  techniques:  predictor-corrector  methods 
which  use  Newton  iteration  in  the  correction  process.  For  a  system  of  m-equations  in 
m-unknowns  such  as  H(y)  =  0,  the  Newton  method  can  be  written  as 

y(s+l)  =  y(s)  .  J-I(y(s))  H(y(s)),  s  =  0,  1,  2  . .  . 
where  (s)  is  the  iteration  number  and  J'1  is  the  inverse  of  the  Jacobian  matrix  9H/3y.  If  this 
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method  is  applied  to  equation  (5.1),  we  obtain  the  sequence  of  approximations 


(s+1)  (s) 

yn+i  -  yn+i 


-  [  I-Atp 


3F(ySi) 

dy 


-1 


]  (jfi+l  "  AtpF(y„+1)  -  q  s  =  0,1,2 


,(s) 


where  I  is  the  identity  matrix. 

Proposed  integration  techniques  to  be  used  in  conjunction  with  Newton  iteration  include  the 
implicit  methods  mentioned  above  and  stiffy  stable  methods  discussed  by  Gear  [41].  These 
methods  are  presented  in  the  predictor-corrector  format  with  a  p-th  order  predictor  formula  of  the 
form 

y  n+l  =  aiyn  +  •••  +&pyn+l-p  +  T]AtF(yn) 


and  a  corrector 

(s+1)  *  *  *  (s) 

yn+1  =  aiyn  +  •  •  •  +  apyn+l-p  +  AtF(yn+i) 

3fe  jfr 

with  a! ,  rj ,  oq ,  and  rj  constants  depending  on  the  order  p.  This  format  differs  significantly 
from  the  more  conventional  methods  in  that  only  one  rate  term  is  used  with  several  previous 
solution  values  rather  than  several  rate  terms  with  one  previous  solution  value. 

All  of  these  predictor-corrector  methods  also  have  infinite  regions  of  absolute  stability  (for 
linear  dynamical  systems)  and  the  use  of  Newton  iteration  provides  for  better  convergence  limits 
than  standard  direct  iteration  techniques.  A  drawback  to  these  methods,  however,  is  that  they 
generally  require  the  calculation  and  "inversion"  of  the  Jacobian  matrix,  for  each  time  step  or,  at 
least,  periodically  during  the  solution  process.  For  only  mildly  stiff  problems,  such  calculations 
may  be  more  time  consuming  than  simply  using  a  smaller  time  step  size  and  direct  iteration.  For  the 
test  problems  considered  here,  Gear  stiffy  stable  methods  performed  better  than  the  other 
predictor-corrector  methods,  with  numerical  results  being  the  same,  to  within  plotting  accuracy,  as 
those  obtained  using  the  forward  Euler  predictor-trapezoidal  corrector  with  time  step  control. 


6.  Test  Cases 

The  computational  methods  and  integration  techniques  discussed  earlier  were  tested  on  a 
series  of  one-dimensional  problems  to  determine  which  techniques  performed  best  and  some  of 
these  results  are  given  in  Figs.  4  -  24.  These  problems  were  selected  from  results  reported  in  the 
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literature,  so  that  comparisons  of  the  numerical  results  could  be  made  with  published  data.  Also, 
along  these  test  cases  are  problems  which  involve  many  of  the  complex  loading  histories  which 
intemal-state-variable  constitution  equations  are  capable  of  modeling.  While  only  two  sets  of 
constitutive  equations  were  considered,  the  Bodner  and  Partoms  equations  in  both  the  original  and 
deviatoric  form  were  used  in  these  calculations.  The  following  shorthand  notation  is  used: 

BPO  -  Bodner-Partom  equations  in  the  original  form 

BPN  -  Bodem-Partom  equations  in  the  deviatoric  form 

H  -  Harts  equations 

TI  -  titanium 

A1  -  1 100  aluminum 

CU  -  OFHC  copper 

SS  -  304  stainless  steel 


In  all  problems  using  Bodners  equations  the  reference  temperature  is  room  temperature,  while  the 
problems  dealing  with  Harts  equations  were  at  temperatures  250°C  and  400°C  for  the  1100 
aluminum  and  304  stainless  steel,  respectively.  The  other  material  data  may  be  found  in  Refs.  [15  - 
27]. 


Test  Problem  #1.  Constant  Strain  Rate  Tests 

Several  constant  strain  rate  simulations  were  performed  under  the  following  conditions: 


EquType 

Mat  Type 

Strain  Rate 

Total  Strain  % 

BPO/BPN 

TI 

3.2E-3/sec 

2% 

1.6E-  3/sec 

2% 

1.6E-4/sec 

2% 

1.6E-5/sec 

2% 

BPN 

CU 

2.0E- 3/sec 

1% 

2.0E-4/sec 

1% 

2.0E- 5/sec 

1% 

H 

AL/SS 

3.33E-3/sec 

1%  -  2% 

3.33E-4/sec 

1%  -  2% 

3.33E-5/sec 

1%  -  2% 

This  large  number  of  constant  strain  rate  tests  were  performed  due  to  the  large  changes  in  stiffness 
of  the  constitutive  equations  for  different  material  types  and  strain  rates.  Note  the  large  differences 
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in  the  results  for  the  TI  material  using  the  original  and  deviatoric  form  of  Bodners  equations,  see 
Figs.  7-8.  Also  note  the  strain  rate  insensitivity  of  the  copper  and  stainless  steel  specimen.  Figs.  9 
-11,  and  non-hardening  nature  of  the  aluminum  specimen  under  these  conditions,  Fig.  10. 


Test  Problem  #2.  Strain  Rate  Change  Tests 

Several  strain  rate  change  tests  were  simulated  under  the  following  conditions: 


Equ  Type 

Mat  Type 

History  of  Loading 

BPO 

TI 

Strain  at  1.6E-5/sec  to  1%  strain  then 
change  to  3.2E- 3/sec  and  strain  to  2% 

BPO 

TI 

Strain  at  3.2E-3/sec  to  1  %  strain  then 
change  to  1.6E-5/sec  and  strain  to  2% 

H 

SS 

Strain  to  3.33E-3/sec  for  0-30  sec  then 
change  to  3.33E-4/sec  for  30-60  sec. 

Results  for  this  type  of  testing  may  be  seen  in  Figs.  12  -  13.  Generally,  this  type  of  loading 
did  not  lead  to  numerical  difficulties  when  changing  from  a  high  strain  rate  to  lower  rate,  but  caused 
some  problems  when  changing  from  a  low  to  high.  This  was  due  to  the  larger  time  steps  allowd  in 
the  initial  lower  strain  rate  simulations.  Also  note  the  asymptotic  approach  of  the  results  to  the 
constant  loading  situation  for  Bodners  equations,  while  Harts  equations  predict  an  almost  immediate 
jump  in  the  response  to  the  constant  loading  case. 


Test  Problem  #3.  Loading,  Unloading,  and  Reloading  Test 
A  single  example  of  loading  at  one  strain  rate,  unloading  into  the  "elastic  region",  and 
reloading  at  a  different  strain  rate  was  performed  under  the  following  conditions: 

Equ  Type  Mat  Type  History  of  Loading 

BPO  TI  Strain  at  1.65E-5/sec  up  to  1%  strain, 

unload  into  the  "elastic  region"  and  then 
reload  at  3.2E- 3/sec  up  to  2%  total  strain. 

Results  for  this  test  problem  using  the  forward  Euler  predictor-trapezoidal  corrector  are 
shown  in  Fig.  14.  The  forward  Euie  predictor  with  time  step  control  from  [12]  experienced 
difficulties  in  the  reloading  part  of  this  problem.  For  this  problem  using  the  old  form  of  the 
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equations  and  the  TT  specimen,  the  predictor-corrector  with  direct  iteration  performed  the  best. 


Test  Problem  #4.  Stress  Relaxation  Tests 

Stress  relaxation  simulations  were  conducted  under  the  following  conditions: 


EquType 

Mat  Type 

History  of  Loading 

BPO 

TI 

Strain  to  2%  at  1.65E- 5/sec  and  hold 
for  70  min. 

H 

AL 

Strain  to  .03%  total  strain  at  3.33E-4/sec 
and  hold  for  10  hours. 

Results  of  this  testing  are  shown  in  Figs.  15  -  16.  For  this  problem  there  were  no  rapid  load 
changes  or  serious  stiffness  problems  and  thus  the  forward  Euler  method  with  time  step  control 
performed  as  well  as  the  other  methods. 


Test  Problem  #5.  Creep  Test. 

Two  creep  simulations  were  performed  under  the  following  conditions: 


Equ  Type 

Mat  Type 

History  of  Loading 

H 

AL 

Rapidly  stress  to  1500  psi  and  then 
hold  constant  over  0-100  hours. 

H 

AL 

Rapidly  stress  to  1000  psi  and  then 
hold  constant  over  0-100  hours. 

Results  are  reported  in  Fig.  17  for  a  l-inch  material  specimen.  This  problem  also  posed  no 
stiffness  difficulties  due  to  the  low  level  of  initial  stressing,  and  thus  the  simpler  methods  arc  more 
efficient  for  this  case. 


Test  Problem  #6.  Stress  Change  Test 

A  single  stress  change  test  was  simulated  under  the  following  conditions: 
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EquType  Mat  Type  History  of  Loading 

H  SS  Rapidly  stress  to  20000  psi  and  hold 

constant  over  0-10  hours,  then  rapidly 
increase  the  load  to  30000  psi  and  hold 
constant  over  10-100  hours. 

Results  of  the  numerical  integration  for  this  problem  are  shown  in  Fig.  18.  This  problem 
having  an  imposed  stress  history  and  relatively  large  hold  times  poses  no  large  stiffness  problem 
and  therefore  the  simpler  methods  are  more  efficient  here  also. 


Test  Problem  #7.  Cyclic  Tension-Compression  Tests 

Several  cyclic  tension-compression  simulations  were  performed  for: 


Equ  Type 

Mat  Type 

Strain  Rate  (=) 

#  Cycles 

±  %  Strain 

BPN 

TI 

3.2E-3/sec 

10 

1% 

BPN 

AL 

2.0E-3/sec 

5 

1% 

BPN 

CU 

2.0E-4/sec 

1 

1% 

H 

SS 

1.0E-3/sec 

5 

0.5% 

(See  Figs.  19  -  22  for  results  of  this  simulation.)  In  this  group  of  problems,  the  constitutive 
equations  were  significantly  stiffer  for  AL,  CU  and  SS  specimens.  Also,  the  region  of  time  over 
which  the  equations  were  stiff  was  large  in  relation  to  the  total  time  and  therefore  the  stiffy  stable 
methods  Gear  perfor.  led  best  in  these  cases.  Note  that  the  number  of  cycles  and/or  percent  strain 
was  in  some  cases  limited  by  the  stiffness  of  the  equations,  so  that  the  predictor-corrector  methods 
using  direct  iteration  were  competitive  in  total  computing  time. 

Test  Problem  #8.  Cyclic  Relaxation 

A  single  cyclic  relaxation  simulation  was  performed  under  the  following  conditions: 

Equ  Type  Mat  Type  History  of  Loading 

BPO  TI  Strain  at  2.5E-3/sec  between  1.25%  total  strain  and 

0.225%  total  strain  over  5  cycles 

Results  of  the  numerical  integration  for  this  problem  are  shown  in  Fig.  23.  Five  cycles  were 
selected  here  because  at  this  point  an  almost  steady  state  was  reached.  For  this  problem  the 
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constitutive  equations  were  only  moderately  stiff  and  the  predictor  corrector  with  direct  iteration  and 
Gears  method,  performed  approximately  the  same. 

Test  Problem  #9.  Cyclic  Creep  Test 

A  single  cyclic  creep  test  was  simulated  under  the  following  conditions: 

Equ  Type  Mat  Type  History  of  Loading 

BPO  TI  Stress  to  325  MPa  at  32.5  MPa/sec  and  then 

to  -225  MPa  at  -32.5  MPa/sec  and  continue 
for  9  cycles. 

Results  for  this  history  of  loading  shown  in  Fig.  24.  This  problem  also  has  only  a 
prescribed  stress  history  and  the  predictor  methods  with  direct  iteration  performed  well  here. 

In  summary,  we  list  for  the  test  problems  considered: 

1 .  The  computational  methods  in  Section  4  produced  identical  results  with  all  the 
time  integration  techniques  discussed  in  Section  5,  with  no  particular  advantage 
of  one  algorithm  over  the  other. 

2.  Gears  integration  package  of  [14]  showed  some  inadequacies  in  strain  control 
situations.  These  difficulties  were  overcome  by  a  slight  modification  of  his 
methods  of  time  step  selection. 

3 .  The  forward  gradient  scheme  presented  in  Section  4  showed  much  the  same  type 
of  deficiencies  as  the  other  forward  type  of  integration  methods.  The  results  for 

a  cyclic  loading  situation  may  be  seen  in  Fig.  6  where  the  time  steps  have  been 

selected  as  a  fraction  of  the  total  strain. 

4.  For  reasons  mentioned  in  Section  5  and  from  results  of  the  test  problems,  no 
method  suggested  in  all  of  the  literature  examined  is  completely  satisfactory.  Of 
those  actually  tested,  the  two  that  performed  best  were  forward  Euler 
predictor- trapezoidal  corrector  with  direction  iteration  and  step  control  of  [12]  in 
conjunction  with  either  the  initial  strain  of  the  initial  strain  method,  and  Gears 
stiffy  stable  methods  wit'  Newton  iteration  corrections  in  conjunction  with  the 
new  algorithm  of  Section  4. 
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7.  Other  Applications 


In  this  final  section,  we  apply  two  of  the  better  techniques  described  earlier  to  some  specific 
engineering  problems  found  in  the  literature.  1  he  methods  used  are: 

1)  the  forward  Euler  predictor-trapezoidal  corrector  with  direct  iteration  and  time  step 
control  of  [23]  in  conjunction  with  the  initial  strain  method  or  the  initial  strain  rate  method,  and 

2)  Gears  stiffy  stable  method  in  conjunction  with  the  new  algorithm  of  Section  4. 


Example  Problem  #1  (Composite  Sheet  using  Bodners  Equations).  This  problem,  taken 
from  [2],  is  essentially  one-dimensional  in  nature.  It  consists  of  a  composite  material  strip  half 
titanium  and  half  copper,  supported  by  two  walls  (see  Fig.  25).  The  specimen  is  rapidly  loaded  at 
the  material  interface  with  a  uniform  pressure  of  150  MPa,  which  is  approximately  six  times  the 
"yield  stress"  of  the  copper,  and  held  constant  over  a  period  of  100  hours.  The  finite  element  model 
and  100  times  the  deformed  configuration  at  time  =  100  hours  is  shown  in  Fig.  26.  The 
one-dimensional  elastic  bar  elements  which  have  been  included  are  essentially  rigid  and  are  used 
here  to  allow  a  homogenous  mode  of  deformation  to  be  modeled.  The  stress  in  the  copper  and 
titanium  components  is  plotted  versus  time  in  Fig.  27.  This  plot  shows  the  rapid  change  in  stress 
from  the  initial  elastic  solution  to  almost  steady  state  values.  Figure  28  shows  the  relative 
displacement-time  curve  for  a  point  on  the  material  interface.  In  this  figure,  Ucj  is  the  elastic 
displacement  of  the  interface.  These  results  compare  well  with  those  in  [2],  even  though  a  different 
constitutive  model,  Bodners  model,  was  used. 

For  this  problem,  the  first  method  of  solution  required  approximately  3  min  20  sec  of  CPU 
time,  while  the  second  method  required  5  min  20  sec  ( on  a  Harris  800 II). 

Example  Problem  #2  (Pressure  Cylinder  usr  s  Equations).  A  hollow  circular 
aluminum  cylinder,  with  an  internal  radius  of  5  inches  and  external  radius  of  10  inches,  is  loaded 
under  plane  strain  conditions  with  an  internal  pressure  of  750  psi.  The  pressure  is  rapidly  applied, 
so  that  an  initial  elastic  stress  distribution  is  present,  and  held  constant  over  100  hours.  The  finite 
element  discretization  for  this  problem  consists  of  five,  eight-node  quadratic  elements  along  the 
cylinder  thickness. 

Plots  of  the  radial  stress,  a:  ,ial  stress,  and  circumferential  stress  distributions  at  various  times 
are  shown  in  Figs.  29-31.  These  results  are  essentially  the  same  as  those  given  in  [5]. 

For  this  problem,  method  1  required  1 1  min  30  sec  of  CPU  time,  while  method  2  required 
22  min  10  sec. 
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Example  Problem  #3  (Perforated  Torsion  Strip  using  Bodners  Equations).  A  rectangular 
perforated  tension  strip  is  loaded  with  a  uniform  tensile  stress  as  is  shown  in  Fig.  32.  The  stress  is 
applied  rapidly,  so  that  an  initial  elastic  stress  distribution  is  present,  and  held  constant  for  a  period 
of  30  seconds. 

For  a  stress  of  100  MPa,  plots  of  the  axial  stress  versus  time  along  section  A-A  are  shown  in 
Fig.  33.  The  advancement  of  the  "plastic  zone"  from  0.5  sec  to  30  sec  is  also  shown  in  Fig.  34, 
with  regions  of  equivalent  nonelastic  strain  at  time  =  30  seconds  shown  in  Fig.  35. 

This  tension  strip  was  also  subjected  to  a  stress  of  125  MPa  and  held  constant  for  30  seconds 
as  above.  The  corresponding  growth  of  the  "plastic  zone"  is  shown  in  Fig.  36.  Comparing  this 
figure  with  Fig.  34,  we  observe  a  large  increase  in  size  of  the  nonelastic  region  for  an  increased 
loading. 

For  loading  up  to  100  MPa,  the  first  integration  method  required  100  minutes  20  seconds  of 
CPU  time,  while  the  second  algorithm  used  49  minutes  20  seconds.  For  the  125  MPa  loading,  the 
first  solution  method  completed  in  226  minutes  10  seconds  and  the  second  method  required  131 
minutes  40  seconds  (on  a  Harris  800 II). 
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Figure  1  Stable  Time  Steps  for  Forward  Euler  Integration  of  Bodner's  Constitutive 
Equations  for  Various  Materials. 


Figure  2  Stable  Time  Steps  for  Forward  Euler  Integration  of  Bodner's  Constitutive 
Equations  for  Titanium. 


Figure  3  Decay  of  Oscillatory  Behavior  for  a  Conditionally  Stable  Time  Integration 

Method. 


Figure  4  Calculated  Stress-Strain  Curve  for  Titanium  with  Time  Steps  Selected 
from  xAe  /  en . 


Figure  5  Calculated  Stress-Strain  Curve  for  Titanium  with  Time  Steps  Se  xted 
from  a  Doubling/Halfing  Technique. 


Figure  6  Calculated  Stress-Strain  Curve  for  Titanium  Obtained  using  the  Forward 
Gradient  Scheme. 


Figure  7  Monotonic  Stress-Strain  Curves  for  Titanium  at  Constant  Strain  Rates 
(Bodner's  Equations  in  the  Original  Form.) 


Figure  8  Monotonic  Stress-Strain  Curves  for  Titanium  at  Constant  Strain  Rates 
(Bodner's  Equations  in  the  Deviatoric  Form.) 


Figure  9  Monotonic  Stress-Strain  Curves  for  OFHC  Copper  at  Constant  Strain 
Rates  (Bodner's  Equations.) 


Figure  10  Monotonic  Stress-Strain  Curves  for  1 100  Aluminum  at  Constant  Strain 
Rates  ( Hart’s  Equations.) 


Figure  1 1  Monotonic  Stress-Strain  Curves  for  304  Stainless  Steel  at  Constant  Strain 
Rates  (Hart's  Equations.) 


Figure  12  Stress-Strain  Curves  for  Titanium  Subjected  to  Rapid  Changes  in  Strain 
Rate  (Bodner's  Equations.) 


Figure  13  Stress-Strain  Curves  for  304  Stainless  Steel  Subjected  to  Rapid  Changes 
in  Strain  Rate  (Hart's  Equations.) 
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Figure  14  Stress-Strain  Curves  of  Titanium  Subjected  to  Unloading  and  Subsequent 
Reloading  at  a  Faster  Rate  (Bodner's  Equations.) 


Figure  15  Stress  Relaxation  Curve  for  Titanium  After  Preloading  to  2%  Strain  at 
1.6E- 5/sec  (Bodner’s  Equations.) 


Figure  16  Stress  Relaxation  Curves  for  1 100  Aluminum  with  Different  Initial 
Hardness  Values  (Hart's  Equations.) 


Figure  17  Creep  Curves  of  1 100  Aluminum  with  Variable  Initial  Hardness  and 
Applied  Stress  (Hart's  Equations.) 


Figure  18  Calculated  Strain-Time  Curve  for  304  Stainless  Steel  Subjected  to  a  Stress 
Change  Test  (Hart's  Equations.) 


Figure  19  Cyclic  Stress-Strain  Curves  for  Titanium  for  ±  1 -Percent  Strain  (Bodner's 
Equations.) 


Figure  20  Cyclic  Stress-Strain  Curves  for  1 100  Aluminum  for  ±  1 -Percent  Strain 
(Bodner's  Equations.) 


Figure  21  Cyclic  Stress-Strain  Curves  for  OFHC  Copper  for  ±  1 -Percent  Strain 
(Bodner's  Equations.) 


Figure  22  Cyclic  Stress-Strain  Curves  for  304  Stainless  Steel  for  ±  0.5-Percent 
Strain  (Hart's  Equations.) 


Figure  23  Computed  Stress-Strain  Curves  for  Titanium  for  Strain  Controlled 

Cycling  with  Positive  Strain  Limits  Showing  Creep  Relaxation  (Bodner's 
Equations.) 


Figure  24  Computed  Stress-Strain  Curves  for  Titanium  for  Stress  Controlled 

Cycling  with  a  Positive  Mean  Stress  Cyclic  Creep  (Bodner's  Equations.) 


Figure  25  Composite  Titanium-Copper  Bar  Loaded  at  the  Material  Interface. 


Figure  26  Original  and  Deformed  Finite  Element  Mesh  for  a  TI/Cu  Composite  Bar. 


Figure  27  Stress  Redistribution  over  Time  for  a  TI/Cu  Composite  Bar. 
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Figure  28  Relative  Displacement  over  Time  of  the  Material  Interface  for  a  TI/Cu  Bar. 


Figure  29  Variation  of  Radial  Stress  in  a  Creeping  Cylinder  in  Plane  Strain. 

Figure  30  Variation  of  Axial  Stress  in  a  Creeping  Cylinder  in  Plane  Strain. 

Figure  31  Variation  of  Circumferential  Stress  in  a  Creeping  Cylinder  in  Plane  Strain. 
Figure  32  Perforated  Tension  Strip  in  Uniaxial  Tension. 

Figure  33  Normal  Stress  Distributions  Along  Section  A-A  in  a  Perforated  Tension 
Strip. 

Figure  34  Time  Dependent  Regions  of  Nonelastic  Deformation  for  a  Perforated 
Tension  Strip  (Applied  Stress  =  100  Mpa). 

Figure  35  Regions  of  Equivalent  Nonelastic  Strain  at  30  Seconds  for  a  Perforated 
Tension  Strip  (Applied  Stress  =  100  Mpa). 

Figure  36  Time  Dependent  Regions  of  Nonelastic  Deformation  for  a  Perforated 
Tension  Strip  (Applied  Stress  =  125  Mpa). 
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Figure  1  Stable  Time  Steps  for  Forward  Euler  Integration  of  Bodner's  Constitutive 
Equations  for  Various  Materials. 
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Figure  6  Calculated  Stress-Strain  Curve  for  Titanium  Obtained  using  the  Forward 
Gradient  Scheme. 
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Figure  7  Monotonic  Stress-Strain  Curves  for  Titanium  at  Constant  Strain  Rates 
(Bodner's  Equations  in  the  Original  Form.) 
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Figure  9  Monotonic  Stress-Strain  Curves  for  OFHC  Copper  at  Constant  Strain 
Rates  (Bodner's  Equations.) 
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Figure  10  Monotonic  Stress-Strain  Curves  for  1 100  Aluminum  at  Constant  Strain 
Rates  ( Hart's  Equations.) 
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Figure  1 1  Monotonic  Stress-Strain  Curves  for  304  Stainless  Steel  at  Constant  Strain 
Rates  (Hart's  Equations.) 
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Figure  1 2  Stress-Strain  Curves  for  Titanium  Subjected  to  Rapid  Changes  in  Strain 
Rate  (Bodner's  Equations.) 
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Figure  13  Stress-Strain  Curves  for  304  Stainless  Steel  Subjected  to  Rapid  Changes 
in  Strain  Rate  (Hart's  Equations.) 


Figure  14  Stress-Strain  Curves  of  Titanium  Subjected  to  Unloading  and  Subsequent 
Reloading  at  a  Faster  Rate  (Bodner's  Equations.) 
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Figure  15  Stress  Relaxation  Curve  for  Titanium  After  Preloading  to  2%  Strain  at 
1.6E-5/sec  (Bodner's  Equations.) 
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Figure  16  Stress  Relaxation  Curves  for  1 100  Aluminum  with  Different  Initial 
Hardness  Values  (Hart's  Equations.) 
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Creep  Curves  of  1 100  Aluminum  with  Variable  Initial  Hardness  and 
Applied  Stress  (Hart's  Equations.) 
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Figure  18  Calculated  Strain-Time  Curve  for  304  Stainless  Steel  Subjected  to  a  Stress 
Change  Test  (Han’s  Equations.) 


Figure  19  Cyclic  Stress-Strain  Curves  for  Titanium  for  ±  1 -Percent  Strain  (Bodner1 
Equations.) 
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Figure  21  Cyclic  Stress-Strain  Curves  for  OFHC  Copper  for  ±  1  -Percent  Strain 
(Bodner's  Equations.) 


Figure  22  Cyclic  Stress-Strain  Curves  for  304  Stainless  Steel  for  ±  0.5-Percent 
Strain  (Hart’s  Equations.) 


Figure  23  Computed  Stress-Strain  Curves  for  Titanium  for  Strain  Controlled 

Cycling  with  Positive  Strain  Limits  Showing  Creep  Relaxation  (Bodner's 
Equations.) 


Figure  24  Computed  Stress-Strain  Curves  for  Titanium  for  Stress  Controlled 

Cycling  with  a  Positive  Mean  Stress  Cyclic  Creep  (Bodner’s  Equations.) 


Figure  25  Composite  Titanium-Copper  Bar  Loaded  at  the  Material  Interface. 
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Figure  26  Original  and  Deformed  Finite  Element  Mesh  for  a  TI/Cu  Composite  Bar. 
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Stress  Redistribution  over  Time  for  a  TI/Cu  Composite  Bar. 
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Figure  3 1  Variation  of  Circumferential  Stress  in  a  Creeping  Cylinder  in  Plane  Strain. 


Normal  Stresses  Along  A-A  (MPa) 


Figure  34  Time  Dependent  Regions  of  Nonelastic  Deformation  for  a  Perforated 
Tension  Strip  (Applied  Stress  =  100  Mpa). 


Figure  35  Regions  of  Equivalent  Nonelastic  Strain  at  30  Seconds  for  a  Perforated 
Tension  Strip  (Applied  Stress  =  100  Mpa). 
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